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ABSTEA.CT 


The problem of the determination of optimal time-invariant output 
feedback controllers for linear dynamic systems "with quadratic cost 
functionals is considered. Tvro distinct cases arise, depending on "whether 
optimization is over a finite or a semi-infinite time interval. 

For the semi-infinite time interval problem a gain initialization 
technique is derived to complement existing optimization techniques . The 
gain initialization technique determines the feedback gains req"uired tq 
(locally) maximize the system stability, A computational algoritlrp for 
the technique is incorporated in a digital computer program, and is uqed 
"(;o stabilize a seven state model of a Satnm V booster rpcket. 

For the finite time interval problem a technique is derived to 
(locally) minimize the expected value of the cost functional. The 
technique uses the concept of Initial State Averaging. A computational 
algorithm is provided and incorporated in a digital computer program. The 
technique is illustrated by three examples. 
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CHAPTER I 
INTRODUCTION 


I.l Motivation 

During the past decade much excellent work has heen done in the field 
of Optimal Control Theory. Many books and a great number of technical 
papers have been published - the Biblipgraphy gf this report cites only a 
small portion of the total. Not all of this work has been purely thepreti- 
cal. A substantial background has been developed in the impor'fcant practice 
area of corapratation and it has for several years been feasible to detersnine 
optimal controllers for a wide variety of physical systoas. 

It would be expected that practical applications of this theory would 
have appeared in abundance.. Optimization is the essence of all ,goo4 
engineering design, and practicing engineers should seize upon any techniques 
which might aid them in their work. Application of the new techniques 
has, however, been disappointingly slow. With a few notable exceptions 
the practical design of control systems has remained based on the transfer- 
function techniques developed prior to the 1960 ’s. 

The explanation of the above paradox is widely recognized. It lies 
in the nature of the word "optimal", which is meaningless without a criterion 
of optimality. By and large, control theorists have used criteria of 
optimality dependent only on the performance of the control system. Design 
engineers, on the other hand, interpret "optimal" as embracing both system 
performance and system cost. .The difference may be best illustrated by a 
simple example. 

Consider the problem of designing a speed control device for a cheap 
movie camera. The objective is to ensure that sixteen frames of film are 





exposed per second regardless of film tension, state of charge of the 
motor "batteries, etc. Ttie could be interpreted as comprising an sample 
of the state-regulator problem in Optimal Control Theory, -with the follp'Wiijig ■ 
method of analysis. The equations governing the film trausportatipp 
determined and linearized to the form 

x(t) = A x(t) + B u(t) (I, Ini') 

■where ^(t) is an WS-vector describing the deviation of the system , 
from its desired s-fcate at time t 
u(t) is an NC-vector defining the control inputs at time -t 
A and B are matrices -whose elements depend on the charaoterr 
istics of the system. 

A cost functional is formulated having "the form 

T 

j = x^(t) p x(t) + [jL^(r) 0 x(‘2') R iiC"?:)] dr (l.l-,2) 

0 

■where T is the d'uration of the scene to be filmed. The objecti-ve is tp 

determine the control input u(t) for t e [ ^ ^ "^bich minimizes -j:he 

cost functional J. The positive semi -definite matrices P and ^ arp 

chosen to penalize any given control system for the deviations it alloif^ 

from the desired state. The positive definite matrix B is a penalty fop 

improvident usage of control power (i.e., battery energy). The solu-j:ipjj 

2 

to the problem defined by (l.l-l) and (1.1-2) can be shown to be 

u(t) = -r"^ D(t) x(t) (I.ln3) 


vhere B(t) satisfies the matrix Riccatl equation 

r»(t) = -D(t) A - A^D(t) + D(t) B^D(t) - Q 





•with 


D(T) = F (1.1-5) 

Defining 

K'^(t) = r“^ B^D(t) (1.1-6) 

and re--writing (I.I-3) as 

u(t) = -K^(t) x(t) (14-7) 

■we see that the optimal control input can he generated as a timervany?.ng 

linear comhination of all of the states of the system. We note in parsing 

that the less-fastidious photographer, mlling to edit a few frames fyom 

the heginning and end of each scene and settle for gopd steady statq 

perfqrmance, would he re-warded hy simplification of the controller? to the 
2 

form 

u(t) = -k"^ x(t) (1.1-8) 

This time- invariant controller would he sauch simpler to mechanize than 
that described by (l.l-Y). 

Mow consider that speed control of the typical cheap movie camera, ie 
achieved by means of a simple flyball governor driving an on/off switch 
between the batteries and the motor. The flyball governor, invented hy 
James Watt in I788, was the first widely used automatic feedback conti'oUp^* 
In its simpler form it suffers from "hunting" about the set ppint. In its 
application to the movie camera it is driven hy the single system state, 
motor speed. The single feedback gain can be considered as constant in 
thb vicinity of the set point. It is clear* that the flyball governor does 
not satisfy the requirements for an optimal controller as defined by (i,1t7) 
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Yet the flyball goyemor is indeed an. optimal controller in the best 
engineering sense. It does an acceptable job at minimum cost. Ihe 
mechanization of the "optimal" controller defined by (I.I- 7 ) on the other 
hand^ would be prohibitively expensive. It would require measurement and 
feedback of all system states. It woirld also require a digital computer 
to calculate and store the time-varying feedback gains for each scene to 
be filmed. 

It shoxild be clear from the above that the allowable degree of 
complexity of a control system is often constrained by consid,erations of 
cost. Weight, reliability and common engineering, sense si.nvi larj.y pften 
dictate sin 5 >licity. What the control system designer wants then is not 
that system which performs in the best possible manner. Rather he wis^^s 
the best possible performance for a given degree of complexity. This 
problem is considerably more difficult than that defined by (l.l-l) and 
(1.1-2). It must however? he solved if the benefits of optimal control 
are ever to be realized for the majority of potential applications. 

This report considers a portion of the general problem of pptimizatipn 
within a given degree of control system qomplexity. Only linear systems 
with quadratic cost funptionals are considered. Such combinations may be 
described by (l.l-l) and (1.1-2). The allowable degree of cpiaplexity of 
the controller is assumed to be time-invariant feedback of the system outputs 
only. The system outputs are those quantities which can be measured? and 
in most cases provide less than a full description of the system state at 
any time. The above limitations are those which are normally applied in 
classical control system analysis techniques. Thus this report seeks to 
optimize the classical controller, without increasing its complexity. 
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T.2 Historical Review 

It is generally conceded that control system design prior to the 

Second World War tras primarily an art. The techniques used were often 

enrpirical and thus confined to particular classes of problenip. This 

situation was largely eliminated by the contributions of Wyquist;^ 

*7 8 

Bode and Evans who established the theory and techniques of control 
system design which predominate in practical work even today. These 
techniques, however, are most suited to the design of relatively siiiple 
systems. They -?fork best for systems having a single input and a single 
output related by a transfer function. Also these techniques aye purely 
analytical and cannot be used to directly synthesize a control system. 

The problem of synthesis was ifirst tackled by Wiener^, who considered 
the optimization of linear filters. This t^ork was extended to control 
systems by Newton, Gould and Kaiser^*^, who, however, still retained the 
transfer-function approach; their technique was to determine the values 
of system parameters required to minimize a cost functional for a specific 
perturbation of the system. It is significant that they made use of the 
Calculus of Variations in their approach to this problem. The extension 

to the classical Calculus of Variations provided by Pontryagin’ s Maximup\ 

11 12 
Principle, and the control problem, framework provided by Kalman finally 

brought Optimal Control Theoiy to maturity. This maturity, combined mth 

the state-variable formulation presented by DeRusso, Roy and Close^ and 

the increasing capability of digital computers, finally allowed the deter- 

inlnation of optimal controllers for a wide variety of systems. 

The early excitment at the disclosure of the ability to compute 

optimal controllers faded when it was realized that such eontropLle^s werq 

very difficult to mechanize. The reasons for this difficulty are twofold. 

First, an optimal controller requires knowledge of the complete state of 
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the system. It may not, however, he feasible to measure all system states. 
Secondly, for finite time duration problems an optimal controller requires 
time-vaiying feedback. 

There are two fundamental approaches to elimination of the first 

- IS. 

difficulty. The approach taken by Kalman and Bucy was to estimate the 

unmeasurable states from the (noise corrupted) system outputs or measurable 

l4 15 

states. This approach was extended by Luenberger and again by Ash . 

Estimation of the unmeasurable states allows the mechanization of a regular 

optimal controller. Note, however, that it addes the complexities of a 

state estimator to those of the optimal controller. 

The second approach to the problem of unmeasurable states is to 

design a controller which uses only the available information, Sych a 

controller will be "sub-optimal" in the mathematical sense of minimizing 

a cost functional such as that described by (1.1-2). It may well, however, 

be optimal in the engineering sense. This approach has been tried by a 

10 

number of people. Newton, Gould and Kaiser essentially advocated such 

controllers, but optimized them for specific system disturbances. (lonsid- 

eration of only a specific disturbance is equivalent to imposing a specific 

initial condition on the system to be controlled. Max-min techniques 

were developed in an attempt to. eliminate the dependence of the solution 

on the chosen initial condition.^^^^'^ In the m^x-min procedure the 

TnaT fTwuTTi system cost with respect to a set of ;initial conc^tions is minimised 

with respect to the system feedback .gains. It still remains to choose an 

l8 

appropriate set of initial conditions. Rekasius circumvented this 

problem by considering the initial condition giying the maximum ratio of 

output feedback cost to allstate feedback cost but was still limited to 

19 

design for this specific initial condition. Levine determined the 
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controller which was optimal in an average sense and computed time varying 

feedhaci; gains for the case of the finite time-interval problem. For the 

semi-infinite time intenral problem he required his feedback gains to be 

20 

time-invariant. Gag sidy considered the minimization of a modified cost 

functional having the form 

T 

If T* Tjl. T T T 

. i ~ - - - - . . 

A A 

where the matricets Q and W were chosen to ensure that the feedback 
gains for the unavailable states were zero. The resultant feedback gains 
were again time-varying for a finite time interval, time-invaxiant for th^ 
semi-infinite time interval. 

It is considered that the approaches used by Levine and Cassidy are 
the most satisfactory of those considered, since they are both conpletely 
independent of Initial conditions. They both, however, result in time- 
varying gains for the finite time-interval ease. For the semi-infinite 
time interval both give constant feedback gains, but use computational 
algorithms which require initialization by a set of stabilizing feedback 
gains . 

The problem, of determining stabilizing output-feedback gains for 

conplex systems has not, to the author's knowledge, been satisfactorily 

21 

solved. Koenigsberg has considered local stability maxima; but his 

22 23 

computational algorithm is felt to be inefficient. Jameson and Davison^ 
have shown that as many system roots as there are feedbaei^i states can be 
arbitrarily determined, but tbe remaining roots are then unconstrained- 

A solution to the problem of determining optimal constant or pieceyisp 
constant feedback gains for a finite time interval has heqn derived by 
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25 IQ 

Kleinman, Fortmann and Athans. As mth Levine, they considered 
optimality in an average sense, hut confined attention to the allstate 
feedback ease. 

1.3 Scope and Contribution of This Work 

Ehis report considers the determination of optimal time-invari^t 
output feedback control gains for systems ■which can be described by the 
equation 


x(t) = A x(t) + B u(t) (1.3-1) 

■where x(t) is an NS-vector describing the state of the system at time 't 
u(t) is an KC-vector- describing the cpntrol inputs at time t 
A,B are time invariant matrices of appropriate order. 

It is assumed, -without loss of genersility, that the system o-u-fcputj?' 
are the first HF states. If this is not the case then a new set pf q-tjate 
variables may be obtained as follows to satisfy the ass'uo 5 )tion. ^uppqse 
tkat the system described by (1.3-1) has MF independent ou-tputs y(t) 
given by 

y(t) = C x(t) + D u(t) • (I.3r2) 


We shall require that 


u(t) = - K y(t) 


(1.3-3) 


■where K is the time- invariant matrix of feedback gains. Thus 

y(t) = a x(t) ( 1 . 3 -k) 

where I is the Identity matrix of order NF. Assuming |^I rf- DK^j to 


be non-singiilar gives 
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y(t) = C x(t) 


(1.3-5) 


where 


C = fl + DK^l"^ C 


'We note that in most cases D will be zero^ leaving C equal to C. We 
now form, the vector x(t) comprising any (WS ~ WF) of the x(t) state 
variables which are independent of the y(t) variables. We note that 
such variables must exist since x(t) spans NS space while ^(t) only 
spans IJP space. Appending x(t) to y(t) gives 


■(■fc) = 



where Ip is a permutation matrix describing how x(t) was chosen from 
x(t). 

« It follows that 


z(t) = x(t) 

Ip 


C C 

z(t) = A x(t) + B u(t) 

Ip ■ Ip 


ihus;, from (I.3-?) 


C C -X c 

z(t) = A z(t) + B u(t) 
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and z(t) is the desired state-vector having its first KF spates as the 
system outputs. 

Thus the control inputs n(t) for the system described by (l,3-l) are 
to be optimized over the set 

u(t) = - x(t) (1.3-11) 


■where K;, - the time invariant matrix of feedback gains. 


■the form 

NC 



k . . 

\ M 

k; = 

0 

0 


is constrained to 


(I.StIS) 


The cost functional is assxuned to be of the form 


T 

J = x^(T) k x(T) + u(r)]dr- (1,3-13) 

0 

where F and Q, are positive sem i -definite 
E is positive definite 

and optimization is over the time Interval ^0 , T j . 

Two separate cases are investigated. For the pase where the terminal 
time T is infinite it is considered that the techniques presented by 
either Levine or Cassidy are adequate. Both these technique^, however, 
require initialisation by stabilizing feedback gains. In praptioe, indeed 
the author has found that feedback gains giving marginal stability may be 
inadequate, due to nmerical considerations. Chapter II of this report 
therefore presents a computational algorithm designed to find feedback 
gains of the form given by (1.3-12) such that the system stability is driven 
to a local maximum. Note that maximum stability is here defined to mean 
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that the least stable of the system, roots is made as stable as possible. 

i 

The algorithm presented is felt to be superior to that giren by Koenigsberg' 
in that it is computa'bionally faster. 

For the case -where the teminal time T is finite^ an apprp^eh is 

19 

taten similar to that used by Levine , exc^t that the feedback gains are 
required to be time- invariant . The resulting theory is presented in 
Chapter III, together -with a, computational algorithm to aU.ow the 
mination of the optimal gains. 

It is felt that this report will help to fill some of the gaps in the 
existing ability to determine practical feedback controllers Tjhi,eh are 
optimal in the engineering sense. 
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CHAPTER II 

THE GAIN INJTIALIZA.TION PROBLEM 


II. 1 Introduction 

The problem of deteimination of the output feedback gains required 
to stabilize a dynamic system is an interesting one in its own right. In 
this report^ however^ we are concerned with the determination of such gains 

as a pre-requisite to the application of optimization techniques. 

19 20 

Both Levine and Cassidy have considered the determination of 
optimal output feedback controllers for linear systems. Both haye derived 
iterative techniques resulting in time invariant feedback gains when the 
system is optimized over the semi -infinite time interval ^0, coj . Each 
technique, however, requires initialization with a set of stabilizing 
feedback gains in order to ensure convergence of the computational algorithm. 

* 26 

Based on the author’s practical experience with C assid y* s technique 
it appears that an extra requirement on initial system stability may be 
added when the optimization algorithms are mechanized on a digital con5)uter. 
It was found that with low order systems (five states or less) marginal 
stability was sufficient tp ensure convergence of the optimization algorithm. 
When a twenty-one state model of the Saturn V booster rocket was considered, 
however, initialization by feedback gains giving marginal stability was 
insufficient to cause convergence. This was judged to be due to numerical 
truncation resulting from the finite word length in the digital computer. 

The problom was solved by computing a new set of feedback gains giving 
better than marginal stability. Re-initialization of the optimization 
algorithm with the new feedback gains resulted in convergence. It ;is 
considered likely that Levine’s, algorithm would have similar cljaraeteristics . 
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It is clear that practical application of Levine’s and Cassidy's 
techniques requires a method for determination of stabilizing output 
feedback gains. It -would be desirable that the gains give more than 
marginal stability, A method is described in this Chapter for deter- 
mination of such gains, when they exist. 


II -2 Theory of the Gain Initialization Technique 

The problem under consideration is as follows. Givon the system 


x(t) = A x(t) + B u(t) (II. 2-1) 

where x(t) is an NS-vector describing the state of the system at 
time t 

u(t) is an NC-vector describing the control inputs at time t 
A,B are matrices whose elements depend on the characteristics 
of the system 

■?d-th 

u(t) = -K^ x(t) (11,2-2) 


where K is a matrix of feedback gains such that 


x(t) = 


where K is constrained to the form 
NC 



(II. 2-3) 


(II. 2-4) 


how should the feedback gain matrix K be modified, within the constraints 
of (II.2-4) to increase the stability of the. least stable eigenvalue of 
(II. 2-3) 



i4 
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The solution hinges on an escpression given by Fadeev and Fadeeva 

for the sensitivity of an eigenvalue of a matrix to a parameter of the 

matrix. Given any eigenvalue K and corresponding row and coliann eigen- 
T 

vectors v and w of a square matrix D, the sensitivity of the 
eigenvalue to a parameter a of the matrix is given by 


do: 


T ^ D 
- la.- 


T 

V w 


(II. 2-5) 


We note that v and w saftisfy 


D V = >s V 


(II. 2-6) 


Dw = X w 


(II.2-T) 


To apply (II.2-5) to the system descr^.bed by (II.2-3) ^d (11,2-4) 
we note that 




h k 


PI 


= [°i 


(II. 2-8) 


p^^ column 


*fcll 

•«fhere b is the q. column of the matrix B. Thus if k . is any 


-a 


eigenvalue of (II. 2-3) and if v^ and m are corresponding row 


and 


column eigenvectors we see that 

th 


column 


I X. 

i_ 

TF“ 


v.'^ f 0 i b ! ol w. 
-XL ? -q. ! J 


pq. 


T 

V. w. 
—1 —1 


(II. 2-9) 
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The sensitivity of -the real part of to the real feedback gain 

k is simply the real part of the ri^t hand jside of ( II. 2-9). 

ir M. 

Suppose now that it is desired to decrease the real part of 
The direction of the steepest descent in feedback gain space is defined 
by the gradient matrix G Tidiere 





( 11 , 2 - 10 ) 


p = 1, IjF 

q — 1^ * * m ) NC 

II. 3 The Gain Initialization Algorithm 

A gain initialization algorithm ha,s been derived from t]ae expressions 
given in ( II. 2-9) and ( II, 2-10), ^d has been incorporated in the digital 
computer program GEADGN. A listing of the program GEADGU is given in 
Appendix I. A description pf the algorithm/program follows. 

The algorithm is Iterative. It increases the stability pf the least 
stable eigenyalue of the system at each step. If, during a giTen interation 
a new eigenvalue should become less stable than the one being operate^ on, 
this fact Is taken into account in subsequent iterations. 

Consider the rearrangement of the variable elements of the feedback 

on 

gain matrix K into an UF x NO vector K . The rearrangement is 
performed by column ordering of K so that 


Sq " ^(q-l)!®’ + p 

The gradient matrix G is similarly rearranged giving 


(II. 3-1) 



^(q-l)NP + p 


( 1 ^. 3 -?) 
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The steepest descent method requires that variations in the gain 

DO DO 

vector K must lie along the negative of the gradient vector G- . 
Thus 


Ox T.® 0^0 

= - p G 


(II. 3-3) 


•where |j, is a positive constant determined by the step sise to be taken 

and is a variation in the gain vector. Suppose now that we ■wish 

no 0 0 

to determine the variation AK in the gain vector K' required to 
cause a small negative change A he{ iii the real part of the eigen- 

value The variation AK may be computed from the relationship 


0^0 T a^rrO 




(11.3-4) 


When the variation "AK lies in the direction of steepest descent 
then^ from (II. 3-3) and (11.3-4) 




■0 -0 


G ARefh.^ 


(II. 3-?) 


whence 


and so 


-ARefx.J 
^ = Mo^mi 2 


H 0 


II 


- ARe { a a 


0^0 2 


(11.3-6) 


(II. 3-7) 


When the variation A Re { is sufficiently small the syst^ 
corresponding to the feedback gains (*^K° + °AK ) jshould have an eig 

value whose real part approximates Re { k + A R®{k^^ . 
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It remains to determine suitable values for the step size A Re { X^.} . 

It is assumed that in most pases a user of GRADGN ’will have liti;le idea 
of what comprises a suitable step size. The step size is therefore vari.ed 
adaptively based on experience ’with previous iterations. The initial step 
size may either be assigned by the user or given a default value of -0.1. 

In either case the initial step size is used as a program criterion for 
termination due to diminishing returns. Should two successive iterations 
of the progr^ fail to inprove the stability by the initial ptep size, 
the program is terminated. It sho-uld bp noted that apart from its pse as 
a terminatipn control the program is quite insepsitiye to the Initial s^tp 
size, due to the rapidity of the step size adaptation. 

The adaptation control for a given iterative step is descril^ed by 
Figure II. 3-1. One of two basic step size a<yustments is utilized, depending 
on the convexity of the trajectory of the least stable eigenvalue. If the 
trajectory is sufficiently convex then a quadratic curve is fj,1:ted to two 
points on the trajectory and the trajectory gradient at one of the points. 

This procedure is illustrated by Figirre 11,3-2. Ihe gain vector corresponding 
to the minimum point of the quadratic curve. i,s cong)Uted and tested fpr 
stability. If the trajectory is insufficiently convex then the step size 
is doubled. One or both of these procedures may be used repeatedly during 
a single iteration. Note that only one gradient vector is computed pqr 
iteration. 

The initiation point for the next iteration is always the most stable 
gain vector found. The starting step size for the next iteration is set 
to one half of the current inprovement in the least stable eigenvalue. 

This aU.ows the step size to be changed by several orders of the magnitude 
in either direction during a single Iteration, while alloT/ing for the 
decreasing rate of Improvement as a stability maximum is approached. 
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Ki^ AK^ - from previous iteration 



= feedback gain vector at coH^iutation point i 

X. = real part of least stable eigenvalue corresponding 
to feedbaqk gains 

= gradient vector at computation point 1 

AX. = step size from computation point 1 to computation 
point i 

FIGURE II. 3-1 Step Size Adaptation Control for the 
Gain Initialization Program GRAPGH' 


















^9 



AX. = step size from computation point 1 to 
computation point i 

X . = real part of least stable eigenvalue at 

congjutation point i 


Quadratic curve fit to variation of X -with A X 
gives 


and 


?\. = a(AX) +PAX+C 


c = X 
b = 


1 

^ X 


i4X 


= 1 


a = 


Xg-b^Xg-c 


Curve has minimum poipt at A X^ -where jgX 

(AXJ^ 


^ ■ 2a 


2(Xg-X^-AX2) 


= 0 


FIGURE 11,3-2 Step Size Adaptation in 

GEADGN by Quadratic Curve Fit 
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!Eie computer program GRADGH incoxporates an option for the degree 
of stability required for termination. Setting the input parameter ISTOP 
to zero instructs the program to maximize stability. Termination is then 
initiated by reaching a point of diminishing returns, as mentioned above. 
If ISTOP is set of 1 a minimum desired degree of stability, STOPR, may be 
input to the program. The program then terminates when the real part of 
the least stable eigenvalue becomes less than STOPR. 


II. 4 Example of the Use of the Gain Initialization Program GRARGN 

The use of the gain initialization program GRADGR is illustrated in 
this Section by output feedback stabilization of a seven state model of 
the Saturn V booster rocket. The model represents the rocket dynamics at 
a point occurlng 80 seconds after lift-off. 

The seven states considered are 


= measured pitch attitude angle 0^^ 
Xg = measured pitch rate 0^ 


R 


^4 


aerodynamic angle of attack a 
first bending mode deflection 7^^ 
first bending mode deflection rate 7^^ 
engine gimbal angle S 

0 

engine gimbal angle rate P 


It is assumed that only the first two states and x^ are to 

he measured. Thus only these two states are available for feedback. 
Control is achieved via an actuator driving the engine gimbal angle at a 


rate proportional to the actuator input. 
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The model may he represented hy 


x(t) = A x(t) + b u(t) 

with 


(II.J+-1) 


:(t) = 


- ^2 





(11.4-2) 

0 . 

1 . 

0 . 

0 . 

0 . 

0 , 

0 .' 


0 . 

0 . 

.2030 

-.6535 

-.0020 

2.568 

0, 


0137 

1. 

.0407 

.0002 

-.0146 

-.0334 

0. 


0. 

0. 

0. 

0. 

1. 

0. 

0. 

(11.4-3) 

0. 

0. 

0. 

- 44.67 

-.1337 

254.6 

0. 


0. 

0. 

0, 

0. 

Q. 

0. 

1. 


0. 

0. 

0. 

0. 

0. 

- 50 . 

-10. 




0 . 


b = 


0 . 

0 . 


0 . 

0 . 


1 . 
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The derivation of this model is described by Cassidy. 

The program GRAD&N was used to determine the feedback gains k^ and 
corresponding to a local stability maximum of the above system. The 
program output is given in Figure II. 4-1. 
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Figure II. 4-1 


Comtputer Prlnt-Out For 
Example Frot)lem. 
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GA'lfl .SEARCH PROCEEDS ALOMG A NEW- GRADIENT,- rTERATlON piUHBER 


-E.IGV,eC' ERRQR messages 
SH l=D.834SE-06 -• 


.ETfeR= 


DIE^0.3576B-06 


EIGENVECTORS CORRE'S TO MRP ' EIGENV'AtUEv 


■ROW REAL .PART 
-0;i4U699E-01 
0,100000DE 01 
• :a/A?.75&P2E-.-Q0 
0.2059637E-03 
f-O.I463000E-01 
- Q,2?»3? , ? 8E;Q9 
•-Q.217140OE-O1 


ROW I HAG PART - 
O.OOOOGOOE 00- 
O.OOOOOODE 00 
Q.OOOQqOOE.OO--; 
O.OCOOOOOE OO 
-O.OQOOuOOE 00- 

O tQcooooo e O P 

o-OQCooooeoo 


.D. lOOOOElOE- Oh - 
0.4234034E 00 . 

-0.2525228E-26 
■ 0,882TS6SE-26'- 

r-Q.l ?;^960 ier-27 

-0;47a7256e-27 . 


• - O'.OOOOOOOE dOO'. 
O.OOOOOOOE DO . 
-Ti-nonnnnnF ^00’ 

, o.oooooooE 00 :• 
■' t),doooooor 00 

' nVnonnoftnF nn 

.--o;oooooooe so ’. 


GRADIENT MATRIX 
0.2729351E-01 


0'.11556l-7E-0l 


STABIUTY INCREASE STEP SIZE ? 


0.-04000060> 


GAINHATRIX K 
_-0.12:42759E 


-0.708895E-01 
-0. 7088956-01 
-0V499273E 01 
-0,4992736. 01 
-0,452182E OD 
-0,3^09.97^00- 
0,2398376-01 


0,6668106.01 
-:0.'6668l06 -01 

_ . 0,’5 0 1'Pl ie , 0 1- ! 

;t-0,S01011£ 0X‘ 

' O.OOOOOOE 00" 

. .0.0000006 on 

0,0000006 DO ' 


STEP S'l.-Zfe IS 00U6L6D 


GAIN HATRIX.-K 
-0.24855186 01 


-0.10523776 01 


I RdQT-S'-r-?’ -'REAt;' PART 
. 1. 3;' 747 3 66 Er:0 l ’ 
-0. 7473666-01 
-0.4985416 01 
-0.49854XE 01- 


•' ' I HAG, PART.- 
■ 0.66530D6 01 
-0, 6653006 ‘Ol- 
0.5020386 01 




FIGURE II. 4-1 (Cont^d) 
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STcP SIZE IS DDUatEO 




-0,^970828 01 
-0.36^»258g 00 
O.I^asllE 00 


-0.5041I7E 01 

o.ocooooe 00 

0»256780E-QI- 


STEP SIZE IS DOUBLED 




ROOTS REAL PART 

-0.9795488-01 

-0.979548E-01,-' 
-0.494167E 01'- 
-0.494167E 01. '■ 
Q.625687E-01 
0.625687B-01 
-a.220044B OO 


STEP SIZE IS DPI 


I NAG. PART 
0.655B47E 01 


0.243981E 00 
-0.243981E 00 


GAIN MATRIX K 
-D.1988414E 


• - v'"- t/' * * ^ 


ROOTS REAL PART 
-3.126488E 00 
-0. 1264885 00 


iMAG. PART- , \‘r" - '.j 
0.642529E oir:^^ 

-0.642529F -'.V? ■ ■■ v~ > 


-0.315896B-01 

-0.907448E-01 


-0.5335046 OD' 
O.OOOOOOE 00. 


sTEP.'‘s-i'2e'-rs dqubl'ed' 






FIGURE II. 4-1 (Cont'd) 
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[ ^GAtN' SEAR cIh. PROCEEDS ALONG A NE« GRADIENT^ ITERAfsION NUMBER 2 


^..'eiGVEC -ERROR -HESSilGES'-”-’ -. ; . ; v - -- > -■>- 

: Si?LsO,a207E-b2 - . rTER= 5 , D1F=.O.OOOOE 

- 4 

00' ^ %) 

1 EIGENVECTORS; CORRES 
U'"' ROW'.REAL PAR'Td ’-. 

TO HRP EIGENVALUE ' ' ' - . 

■- ■ROU '-TMAG ■'PART..'’: ' rtfU RFAt' PART- 

: . . ' 

cm-- rHAG p-ART-"; -' *• 

p' 0,1000000E 01, 
1 , 0.1029001E 00 

I : -0,-974A528E 00 

'O.-OOOOUOOE QO • 0.2210334E' 00 
O.OCOOOOOE OO ■ 5-0, 1373 143 E-01 
o.flooooooe 00 n-ioononop or 

o.ooooaoOE 00 ' ' 

O.OOOOOOOF 00.' 

o.ooooaooF o'n^ ■ - 

.-rOV1416449E--01 
f ' -0, l’5E9490E-Cf2 • 

-b.l&85rilE-01 

•O.OOOOOOQE 00, ■ O.-OSSOIOSE- 00= 

i'.iO.OOflOdOOE-OO - r-0'i6l37960E-01 

■b.-OOGOODOE 00 • .0.1733184E 00 

. -OiOOOOOOOE 00 ■ . .; 

= O'.OOOOOODE- 00---. ; 
O.QOOOOOQE: 00 ‘ ^ 

. -0.1896896E-02 - 

i 

O.OQOOUOOE 00 -0,10767aie-01 

O.OOOOOOOE DO i 

1 ' • j t 

.'..GRADIENT MATRIX . 
*- ^^--0..543,U92E-03 

0. 3374061 E-»0'4 - 

i ... 



STEP.; SIZE- ?, 0.2,42 754 60:.. ■ 


s' ‘ 

. GAT'hl. MATRIX-- K 
V 0.40547-7 3 B 03 

-0.4449835E.-07- , - 



, r. ,-V V ' c*' ’* * *- 


FROOTS real flART 

IHAG. PART 


i ;^^t)..58&954a-0l 

1 ' -:;,~0.586,954E-. 01 

■ - - ' 0 . 13 7404E -.0 1 ■ 

0.6791D2E 01 . 

‘-0.679102E01 

' '0.-644-78 4 E 01 • . : 


\ .0V1374O4E 01 

1 -• --0.3133'56E 01' 

f 0. 19943, 7S 01 

V0.644784E 01 
O.OOOOOOE 00 ' 
O-OOOOOOE 00 

.. V 

f.) . .,,_--O.,4'4^96'4E-0p 

. .rO.pjOOOOOE 'op ■- ; 

■ ■ • , 

CONVEX' CURVE FIT TO 

STABILITY LOCUS 

-w'.' - 

GAIN'-HAJRIX .'.K' 
.rV0VV620357E;O2 . ' 

Tt0.-I829822E D2 ■ " . ' . . ' . 

'- - - -:u 

■ - '-i 

ROO, TS A "ii'M'iilf R'AR'J. 
i'V -V-^o.483643H''01^- 

V 'Vi,: .IMAGi -PARTV ' 

- ■ ■ 0.54991'7E' 01 " 

' V'\: 6'f 

■ - . -0.483643S 01 

-0.307465E-01 
-0. 307465E-0r 

-0.-5'4991T£: pi, 

0.61^~206E- 01 

-0.612206B, 01 . - . 

i 

- -L 

^:----0>.l?7,808E-;iOO 

-OiasTsOae -oo’ 

i • -O. I24568E OO 

-•-••.;-oi-395643E :ool ■ - - --v .: 

•' '■ -0.395643 E 00 - 
• O.OOOOOOF Oil 




.'I'. V 


FIGUBE II. 4-1 (Cont'd) 
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FIGURE II. ^-1 (Cont'd) 
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FIGURE II « 4-1 (Cont'd) 
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In order to verify that GEADGN did indeed attain a stability maximum, 
the eigenvalues of the system described by (ll.i<— 1) to (II. 4-4) were 
obtained for a matrix of feedback gains. Ihe values of the real parts of 
the least stable eigenvalues were cross-plotted against the feedback gains, 
resulting in the egui-stability contours shown in Figure 11.4-2, Ihe gain 
trajectory produced by GRADGW is superin^iosed on Figure II. 4-2. The 
steepest descent nature of GRADGN, and the fact that a stability maximum 
was achieved {within the limits of the diminishing returns termination) 
are evident from Figure I I. 4- 2. 

11,5 Comments on Gain Initialization 

The digital computer program GRAIXiN clearly provides a way to improve 
the stability of a system by output feedback. As anch it is a useful 
adjunct to the output feedback optimization techniques of Cassidy and 
Levine. It must be remembered, however, that GEADGN is designed to find 
purely local stability maxima. Thus it may on occasion fail to find 
feedback gains giving sufficient stability for initialization of the 
optimization algorithms, even though such gains actually exist. In an 
attempt to ameliorate this problem, provision has been included in GEADGN 
for initializing its gain search at any desired location. Thus any desired 
volume of gain space may be searched by initializing GEADGN at a suitable 
number of discrete points, 

Since GEADGN is designed to perform a similar function to Koenigsberg? s 
algorithm, a coDparison of the two technique's may be pertinent. Both are 
gradient techniques and both use adaptive step size variatipn, although 
Koenigsberg's adaptation procedure is much sicpler than that used in GEADGN. 
A major difference occurs in the eonputation of the gradient. Koenigsberg's 
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gradient confutation is "based on Reddy's expression 


^ Tr (conjoint ^Aj 

^ a ~ Tr (conjoint 


^A 


]) 


W) 


(II. 5-1) 


X = X, 


for the sensitivity of the eigenv^ue X^^ of the matrix A to the matrix 


parameter o:^. The conjoint is defined b 3 '‘ 


conjoint j^Aj = adjoint j^A - X I j 


( 11 . 5 - 2 ) 


The egoivalent expression used by GRADGN is (from ( II. 2-5)) 
T ^A 




- 


w 


^ a 


P 


T 

V w 


(II. 5-3) 


■where v and w respectively are the row and column eigenvectors of the 
matrix A corresponding to the eigenvalue X^^. Equation (II. 5-3) ma3'‘ be 
expanded to 


US US 



I 

i=l 

E 

j=i 

V. 

1 


II 

|/0 

US 

I 

us 

I 

V. 

1 

w. 

3 


i=l 

j=l 



which may be reordered to 





US 

NS 




I 

j=l. 

E 

i=l 

w. 

3 

V [-^1 

^Op 

US' 

I 

US 

E 

w. 

3 

V. 

i 


j=l 

i=l 




(11.5-.^) 


(II. 5-5) 



3 ^ 


Inspecijjxon of (H.5-5) shows it to be the same as 



( 11 . 5 - 6 ) 


Now the determination of a colvonn eigenvector of a matrix may be 
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achieved by the computation of only one column of the conjoint. Van Ness' 
EICxVEC^ used for eigenvector coD 5 >utation in G-RADGN, give^ both row and 
column eigenvectors with about the same computational effort as is reguifqd 
for only a column eigenvector. Thus (11.5-6) is computationally more 
efficient than (II. 5-1)- However (11.5-6) is sickly (II. 5-3) with each 
vector inner product replaced by the trace of a vector outer product^ so 
that (II. 5-3) is clearly preferable to (jI.5-6) especially for systems of 
hi^ order. 
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CHASTER III 

THE FINITE TIME INTERVAL PROBLEM 


III.l Introduction 

In this Chapter we are concerned with the determination of the 
optimal time invariant output feedback controller for the system 

x(t) = A x(t) + B u(t) (Ill.lrl) 

where x(t) is a NS vector 
u(t) is a NC vector 
with the cost functional 

T 

J = x"^(t) F x(t) Q x('Jr) + u^(r) E dT: (lII.lj-2) 

0 

where T is a fixed finite time and F, Q, R are suitably positive 
(semi) definite. Thus we requi3?e 

u(t) = x(t). (111,1-3) 


and constrain K to be of the form 


NC 


K = 





(III. 1-4) 


Two itens set this problem apart from the norma^ linear- quadratic 
state regulator problem. These are the requirement for output feedback 
control and the requirement for the time-invariant feedback gains, 14 is 
remarkable that each of these requirements has been individually satisfied 

19 

by different investigators, but by the use of similar techniques. Levine 
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used the Initial State Averaging (I.S.A.) technique to determine the 

optimal time-variahle output feedback controller for the system (ill. 1-1) 

with the cost functional (ill. 1-2). On the other hand Kleinman, Fqrtmann 
25 

and Athans used I.S.A. to determine the optimal time-invarianl allstate 
feedback controller. Ihis report considers the use of I.S.A. to solvq 
both requirements simultaneously. 

The theory imderlying the I.S.A. approach Is presented in Section 

III. 2 

III. 2 Theory of the Initial State Averaging Technique 

The solution to the problem posed by (III.l-l) to (Ill.lrS) is 
underdefined in the sense that the optimal controller is a function of the 
undefined initial condition x(0) . In general each different initial 
condition requires a different set of feedback gains. Since no single 
controller can minimize the cost functional for all initial conditions 
it seems lifise to seek that controller which minimizes the e:^ected value 
of the cost functional. This is the rationale behind I.S.A. 

From (III.ItI) and (llI.l-<3) we see that 

x(t) = A x(t) (lll.2rl) 

■?diere 

A = [ A - Bk"^] (111,2-1) 

Thus for a given initial condition x(o) 

x(t) = j2S(t) x(0) (111.2.3) 

where 0(t) is the state transition matrix of 5 and is flefined by 
0(t) = A j2S(-t) 


(iII,2t4) 
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0(0) = X (III. 2-5) 

Since we are dealing only with linear time-invariant syjstems we have^ 

0(t) = exp l^^t] (111.2-6) 


Eepeating (ill. 1-2) 

T 

J = x^(T) F x(T) + ^ + ^('t) R u(r) d?; (III. 2-7) 

0 

Using (III. 1-3) reduces this to 


T 

J = x^(t) F x(T) + J x^{r) I^Q + KPk'^J x(r) dt 


(III. 2-8) 


whence, by (III. 2-3) 


T 

= x(0)^ 0*^(1) F 0(T) x(0) + J k'^(o) 0^('T) 


0(r) x(0) dt 


(III, 2-9) 


Moving the time invariant x(0) outside the integration reduces (ill. 2-9) 
to 

T 

J = x(0)^ 1^ 0^(T) F 0(T) -f- J 0“(r) [q -i- ICRK^]0(r) d^jx(o) (III. 2 - 10 ) 


0 


which may be rewritten as 


J = X (0^) F x(0) 


(III. 2-11) 
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■where 


T 

P = 0^(T) F 0(T) + J [q+ KRK^] 0(r) d'2? 


e[j] = E [x'^(O) P X (O)] 


- -NB N3 

i/\ ' ^ I Z ha 

Li=i 3=1 


ITo'w consider the x^(0)^ i = 1?3 to be random -variable whose 

distributions are dependen,t on the distribution of the initial conditions 
in state space. Then Ejjj is given by (lll.2-l4) as the expected value 
of a finite sum of random variables. But the expected v^ue of a finite 
sum of random variables is equal to the sum of the- expected values . Thus 

US US 

Z Z ho h<°>] 

i=l 3=1 

Rearranging (ill. 2-15) gives 


US US 


iA- Z Z ®[h3h<°)h"”] 

i=l 3=1 


It is clear from (ill .2-12) that the P. . in (lll.2-l6) |io not depend on 

the x.(0) or the x.(0). Thus (lll.2-l6) can be rewritten as 
^ 3 

NS NS 


Z Z ho®h<p’ht°)] 


i=l 3=1 
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Suppose now that each of the x^(0) has zero mean and that we can estimate 
the covariance matrix V of the initial conditions in state space, i.e., 




(111,2-18) 


Then, from (III.2-I7) and (lU.2-18) 



Z 


NS 


I 


P. . V. . 

ID ID 


(III. 2-19) 


whence 

e[j] = Tr [pv] (III. 2-20) 

Now consider that the sj^-stem described by (ijl.l-l), (ill. 1-3) and 

(III. 1-4) is linear. Thus superposition holds and the system response can 

he scaled -up or down with the magnitude of the Initial conditions. This 

implies that the response of the system to an initial condition lying on 

the surface of the unit hypersphere in state space characterizes the 

response to any colinear initial condition state vector. Furthermore it 

is clear from (III. 1-2) that if any initial condition is ratioed hy an 

2 

amount r then the cost is ratioed by r . These two facts imply that 
the costs associated with the set of initial conditions lying on the surface 
of the unit hypersphere is state space can he used to conveniently charac- 
terize the costs for all initial conditions. It is shown helow that with 
little loss of generality attention may be confined to a uniform distribution 
of the initial conditioi;LS on the surface of the unit hypersphere. 

Suppose that the covariance matrix V is non-singular mth inverse 

~1 -1 

V . Suppose also that the square root of V exists, and let it be 

denoted by Consider the description of the system, (ill. 2-1) in 

terms of a new set of state variables y(t) defined by 



4o 


y(t) = x(t) 


Then (ill. 2-1) 1)60031163 


y(t) = A(NS)^/^ y(t) 


i.e.. 


y(t) = V'^/^ 1 y(t) 


which is similar in form to (lll.2,-l). 

Also we find that 

E[y(0)/(0)] = e[ x(0) ^"^(0) 


Thus 


P NS 


E[y(o) /(o)] 


(NS)"^ E; 


XJ 


NS 


I I [ 

-p=l q=l 


V 




iP 


X (0) 

P 


X 




V 


- 1/2 


Interchanging finite snin and esqiectation gives 


NS NS r 


• 1 ±N»- 

Ey(0)/(0)J^ =-(NS)"^ ^ 

^ XJ - 


E 


P=1 q=l 


x ^( 0 ) x ^( 0 ) 



NS NS 


p=l q.=l 


(III. 2-21) 

(III. 2- 22) 
(III. 2- 23) 

(111.2-^4) 

" 

(111,2^25) 


(111.2-26) 

,(111.2-27) 



lj.X 



V j 

id 


(III. 2-28) 



(III. 2 - 29 ) 


where I is the identity matrix of order US, 

Thus the covariance matrix of the initial state vector y(0) XP 
compatible with a uniform distribution of the y(0) vector on the surface 
of the unit hypersphere. Since such a distribution is convenient and since 
we can, within the limitations of the non-stringent assumptions made above, 
transform our initial state vector so that it is compatible with this 
(iLstrihution, we shall hereafter assume that 



(III. 2 - 30 ) 


Thus, from (ill. 2-20) we see that 



(III. 2 - 31 ) 


and our I.S.A. optimal control is that which mirdmizes 
P is given by (ill. 2-12). 



where 


III . 3 The Initial State Averagipg Algorithm 

The objective of the I.S.A. algorithm is to determine that gain 
matrix K ^ich minimizes the average cost of control as given by (III. 2-31) 
under the constraint that K must be of the form given by (III. 1-4). 

Thus we wish to minimize the gain function GF given hy 

T 

GF = ^ tr r 0^(T)F 0(T) + J 0^(r) + KKK^ ] 0(T) dr 

^ n 


(III. 3-1) 



Now for any two matrices A and B of suitable dimensions we hare 


the trace identity 


tr I^Ab] = tr Bij 


(Ill,3i2) 


Using (ill. 3-2) and the fact that the integral of a matrix is the matrix 
of the integrals of the elements reduces (ill. 3-1) to 


GF = tr F 0(T) 0^(T) + (Q + KEK^) 


T 

Jp(r)f(r)dzj (III. 3-3) 


Since the variable elOT.ents of the ga^n matrix K are unconstrained the 
necessary conditions for K to minimize GF ar^ given by 





0 


(III. 3-^) 


i = 1, . . . ^ NF 
d = 1, .... WC 


Expanding (ill. 3-4) gives the necessary conditions as 


tr 


F 


5>0 r^K 


'T' ^ K 

R K + KR - 


'■i-d 


T 


41^1 



(in.3-^) 


where the arguments of the 0 functions have been dropped for clarity. 
Now for any matrices A and B 



and 


tr 


[a+b] = tr [a] + tr [b] 

Thus in (III. 3 - 5 ) 


tr 


\ 1 ±. 

ik. . 

L ID 


0^+0 


^0 




h k. 




= tr 


F 




F 0 


^0^ 

^k. . 


hy (III. 3 - 7 ) 


= tr 




by ( 111 , 3 - 6 ) 


tr I F 0 ^’ + tr F -||— 0 ‘ 


A0_ 




10 


by (III. 3 - 2 ) 


= 2 tr F 0"^ 


10 


Similarly 


tr 


^k. . 

10 


Ek"^ + KR ^ ^ 


^k. . 

10 


“I T 

1 / 


0 0^ dr 


= 2 tr 


j 

10 ^ 


h3 

(III. 3-7) 

b] 


(111,3-8) 

0 <f dr 

(III. 3-9) 


and 
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tr 

[ci+kkk’^] f 



= 2 tr 

[q+ Ok’^] j 


L- 0 

O-d 

Id J 


L 5' Id J 


(III. 3-10) 

Tims the necessary conditions of (ill. 3-5) reduce to 


W = 0 

•where the MI’ x HC matrix ¥ has its elements defined hy 


(III. 3-11) 


w. . = tr 
Id 




+ 


T 



+ 


Id 


RK 



(III. 3-12) 


Equations (III. 3 -II) and (111.3-12) represent a set of ME x EC 
simultaneous equations in the ME x MC variable elements of the gain matrix 
K. Solution of these equations gives candidates for the ot^imal controller. 
The solution is performed hy Mewton-Eaphson iteration. 

0 D 

Consider the rearrangement of W into an JiE x WC vector W 
hy column ordering 


i . e . j 


w. . 

Id 


¥ 


MF(a-l) 


+1 


(III. 3-13) 


The variable elements of K may be similarly ordered, giving the gain 

0 D 

vector K •where 


K. . 

Id 


0 a 

K 


ME(d-l)+i 

Suppose that °K(n)° is the gain vector at the iteration of 

the Newton- Raphson procedure^ and. that W(n) is the corresponding 


(HI.3-1^) 





vector of necessary conditions. The Newton-Eaphson technique computes 

a . Ji 

K(n+1) , an, inq)roved estimate of a gain vector satisfying the necessary 
conditions, hy 

°K(n-fl) = °K(h)° -n °W(n)° (ill, 3-15) 

■where ^ (ii^) is the (I'fl' x EC) x (MF x EC) gradient matrix given by 

V7 ^”w(n)“ 

= ^Hx(n)r (111.3-16) 

and n is a convergence factor . 

It can be seen from (111.3-13) a^d (III, 3-1^) that computation of 

y(n) is equivalent to confutation of for g = 1, ..., NF, 

pq 

h = 1, ..., EC; p = 1, ..., EF; q = 1, EC, From (III.3-I2) we 

have 


il 
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0 pq gh 
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^ L pq 
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^k 

0 P‘1 


pq 


(III. 3-17) 
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Clearly 




A k „ ^k , 
m gk 


0 


(111.3-18) 


Also, ky (111.3-6) and (ill, 3-2) 
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(III. 3-19) 


so that (III. 3-17) reduces to 




gk 




- tr 


pq 
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A. 


^k ^k 

pq gk 


s . 2 jT 


^k i k 


pq gk 


+ [q + kre"^] 
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^ k ^k , 
pq gk 
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pq 
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(III. 3-20) 


From (III. 3-3), (III. 3-12) and (ill. 3-20) it can be seen that com- 

0 0 

putation of the cost GF, the necessary conditions vector W and the 
gradient matrix ^ for a given matrix K is dependent on the computation 



of a ntmiber of functions of the state transition matrix 0(t). Basic to 


evaluation of these functions is evaluation of 


^ 0(t) 




and 


PI 


^k , 
pq gk 


for p=l, ...jlTF; q=l, . . g=l^...,MF; h=l;,...^WC. 

Consider the state transition matrix 0^^(t) corresponding to a gain 
matrix K given by 


K - 


K. + 


^k 


pq 


k + 

pq 


^k 


S k 




gk 


(III. 3-21) 


where ^ ^ ^gh small. By (ill. 2-4) 0g(k) can he computed 

as the solution to 


d 0A(t) 
dt 



(jII.3-22) 


with 


0a(o) = I 


(III. 3 - 23 ) 


Thus 


d 0g(t) 
dt 






^k 


gk 



0 |^('t) 


(III. 3-24} 

Treating the second term in (III. 3-24) as an external input and 
using (ill, 2-6) 


0A(t) = exp 


[Fa- 


ll 


«o’']*] - / 


exp 


[[- 


BK, 
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] (t-r)j 


B 






dr 


(111.3-25) 



U8 


We now substitute the function described by (III. 3-25) for 0|^(T) 
in (III, 3-25). This is similar to the method for computation of the matjizant 
given in DeRusso, Roy and Close^^ and results in 


0g(t) = (t) 


0^(t-r)B[4f- Xk + 


pq gk. 


0 0° Lpq.^^ gh®J J 


where (t) is the state transition matrix? corresponding to the gain 
matrix K^. Thus, to second order in ^ ^pq. *^^gh^ 

t 

0 X 0 L-no ^ 


[V">- 


Expansion of (III. 3 -2?) gives 


^(t) = 0 (t) - / 0 (t-T) -0 (r) dr 


^0 ^ 


■pq ^^0 


p T 

- r 0 (t-r) B ll- 0 it) dr 
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t 

u 

I- 0 


® 41 
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pq, ^ 0 pq 0 
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(III. 3-28) 


Fow consider that 0(t) can he expanded in a Taylor's jserie? h^out 


Kq giving 


0g(t) = 0j^ (t) + 


^0(t) 


ak 


0 pq 
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^ 0(t) 
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+ terms of higher order 


(111.3-29) 


Conparison of (111.3-28) and (111.3-29) gives 
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^0(t) 




pq 


and. 


T, 




^0 ^ 


® it- 

pq 0 


(III. 3-30) 


^ ^0(t) 

^ ^ T-. 1 ***/\ 

pq K. 0 ^ 


•JO pq. 5 0 ^0 
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* ' ® 4i^ 

0 ° ^ 


'Z 
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: Cr-=) 41 ^ 0 (s) fis 41^ 

0 P<1 P 


(III. 3-31) 


The result given "by (III. 3-30) agrees -with q, similar exprepsipn 

31 

presented by Levine and Athans for the semi-infinite time interval problem. 
We note from (ill. 3-30) and (ill. 3-31) that 


^^0(t) 

^ k ^ 

pq gli 


.T 


K, 
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pq sh 




(III. 3-32) 


The functions of the state transition matrix to be evaluated are 
now seen to be given by 
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10 
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^k 


pq 


0(T) f fir) -If- 0^(T-r) dr 

n p<i 


(Hi. 3-33) 
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(111.3-34) 
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Ak ak , 
Pq gh 
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(III. 3-35) 




0^CC-s) ds d2r 


(111.3-36) 
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(III.3-3T) 


and 
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Ak ik , 
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^0(s ) ^k ^ 0 ^(s ).AK 
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pq gJi 


^ k 
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(111.3-38) 


Consider (III. 3-33) as Taeing representative of (ill. 3-33) through 
( 111 . 3 - 38 ) • Supposing ^ to have distinct eigenvalues vre can evaluate 
0(t) hy 


0(t) = M exp [At] M"^ (III. 3-39) 

■yrhere M is the modal matrix of eigenvectors of S. 

A 

A is the diagonal matrix of eigenvalues pf A 
Substitution of (ill. 3-39) in (III. 3-33) yields 
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pq 
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(III. 3-40) 
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ITow 




— is simply an NS x NC matrix having a 1.0 in it pq. position 




and zeros elsewhere. Thus the ij element of 0 




m 


is given hy 
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^ k 
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( 111 . 3 -^ 1 ) 


where X 




is the n. 


th 


eigenvalue of A. 


Performing the integration yields 


0 ^0^ 
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3-0 
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^1"’ ^2.^ ^3 


-1 -1 
M. M M M B M M. 

lUi n^n^ HgUg pUg nj^q n^Uj^ 


\ (Xn^+Xn-) (Xn^+Xn^) 

• (^nT+Xn^) '13 '15' 

Iff ^ 1 5 ^ ® -e - 

nc=n_ n_.0n_ Xn_ - -Xn_ 

5 3 5^^ 3 3 5 


( 111 . 3 - 42 ) 


Thus (ill. 3-33) can be evaluated in terms of the eigenvalues and 
eigenvectors of A. Equations (111.3-34) throu^ (ill. 3-38) can be similarly 
evaluated. The terms can then be combined by (III. 3-3), (ill. 3-12) and 
(III. 3-20) to give the average cost GE, the necessary condition^ vector 
and the gradient matrix respectively. 

The only other item required for application of the Newton-Eaphson 
iteration technique described by (III. 3-15) is a value for the conyergence 
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factor [i. Determination of a suitable value of the convergence factor 
is described in Section III. 4. 

III. 4 Mechanization of the Initial State Averaging Algorithm 

The algorithm described above for computation of optimal constant 
output feedback gains has been mechanized in the digital computer program 
ISAFT. The program was written in Fortran. IV for an IBM 360 / 5 O computer. 

A listing is given in Appendix II. The program coiiprises a MAIU and 
several subprograms. The functions of the major subprograms are described 
briefly below. 

MAIN: performs most input-output functions and acts as the controller 

for the other subprograms. The MAIN subprogram also computes a suitable 
value for the Newton-Raphson convergence factor p (see (ill. 3-15) and 
controls program termination. The operation is illustrated by the data 
flowchart given in Figure 111,4-1. 

Con 5 )utation of the convergence factor p is performed adaptively. 

A trial step is made with p set to 1.0. The value of the gain function 
GF for the new gains is compared with its former value. If a decrease 
has occurred the program proceeds to the next iteration. Otherwise the 
value of p is halved. This process is repeated up to five times. If 
no io^jrovement in the gain function has resulted by the fifth cycle, i.e., 
p = 2 , the program is terminated due to poor convergence. 

Assuming no termination due -to poor convergence, normal program 
termination will occur when the feedback gains approach their final values . 
The criterion for this termination Is that the proportional change in 
each individual gain during a single iteration should be less than an 
amoxmt GNSTOP. The parameter GNSTOP is input by the user. 
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Subroutine S 
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I ~ — I 
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p^IGURE Data ?''lowchart for Digital Computer 

Program IS AFT 
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STRAM: computes the eigenvalues and the model matrix of the matrix 

The subprograms VECT, HSBG and ATEIG are used in the eigenvalue 
determination. The subprograms MSQ and EIG\TEC are used to determine the 
model matrix and its inverse. 

FEFN: confutes the functions of the state transition matrix 

described by (ill. 3-33) through (111,3-38). 

GAmrU: computes the value of the gain function GF at the 

beginning of each iteration - see (ill. 3-3) • 

0 0 

WESCON: computes the necessary conditions *vector ¥ - see 

(III. 3-12). 

GRADWT: computes the gradient matrix ^ - see (III. 3-20). . 

INVERT: inverts the gradient matrix 

WEWEIT: perforins one iteration of the Newton-Raphson algorit^ 

see (III. 3 - 15 ). 

GAIN2: computes the value of the gain function GF at the end of 

each iteration for use in the determination of a suitable convergence 
factor p - see (III. 3-15). 

Ill . 5 Examples of the Use of the Optimization Program ISAFT 

Three exaa^iles are given to illustrate the use of the digital 
computer program ISAFT. In keeping -with the expository nature of this 
Chapter the examples are all short and are designed to illustrate different 
aspects of the computations. 

Example 1 

This example considered the case where all system states were 
• available for feedback. The system equations were 
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The control was constrained to he of the form 


u(t) = x^(t) - (III. 5-3) 

It was assumed that there was no prior knowledge about the system^ 
and the initial feedback gains were set to the default value of zero. The 
computer print out is given in Figure III. 5-1. The optimal feedback gains 
were computed to be 


2.30 

(111.5-4) 

3.42 

(III. 5-5) 


To verify that the values given by (ill. 5-^) and (ill. 5-5) were 
indeed optimal the results were checked by independent means, A matrix of 
feedback gains in the vicinity of the solution gains was chosen. At each 
set of gains the system described by (ill, 5-1) and (ill. 5-3) was simulated 
over the optimization time interval^ starting from a number of different 
initial conditions, and the costs were confuted numerically using (ill. 5-2). 

The average cost was then con 5 )uted for each set of feedback gains. The 
initial conditions were chosen to be equally spaced around the unit circle 
in state space. Such an independent check will hereafter be referred to 
as an I.C. Cross-plotting the I.C. data resulted in the equal-cost contours 
shown in Figure III. 5-2, The gain trajectory produced by ISAFT is superimposed 
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on Figure 111,5-2. It is apparent that ISJiFT determined feedback gains 
giving a lainimuiii average cost. 

To provide a con^arison for the results produced by ISAFT, the 

problem defined by (III.5-I) and (ill. 5-2) was solved for the optimal time 

varying feedback gains by the normal procedure of Eicatti matrix backwards 
2 

integration. The resulting feedback gains are compared in Figure III. 5-3. 
with those produced by ISAFT. The average cost for the optimal time 
varying gains was coniputed by numerical integration from a number of initial 
conditions. The average cost for the time-varying gains was I.O3O compared 
with 1,077 for the time invariant gains. 

Example 2 

This example used the same system equations and cost functional as 
Example 1 - see (ill. 5-1) and (ill. 5-2) respectively. The feedback structure 
however, was constrained to feedback of only the first state i.e., 

u(t) = - k^ x^(t) • (111.5-6) 

A rather inaccurate guess at the optimal feedback gain resulted in 
the computer print out given in Figure III. 5-^* The use of the convergence 
factor is evident. 

Figure III. 5-5 compares the results given by ISAFT with those 
determined by an I. C. It is again clear that ISAFT determined a feedback 
gain giving a minimum average cost. 

The minimum average cost for- the single time-invariant feedback was 
1.219 compared with 1.077 for allstate time invariant feedback and I.O30 
for allstate time varying feedback. 

Example 3 


This example was chosen to check the effect of complex eigenvalues 
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rO«lLOJ!CCgp 01 ~C.2179‘i49C 01 


AVERAGE COST = 0.2160.^570.. 01 


NGCESSARV C ONDI TIONS VECTOR 
0.27291080 00 


GRADIENT 'haTr’ix 
0.98928970-02 


INVERSE GRADIENT MATRIX 
. . J<U.082«.9jJ. . 


FIGURE 111.5-4 CConf d) 
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IT6RATIQN NUMBER '. 



SYSTEM, EIGENVALUES. 

. -0;62787600 0,1.,. -O.OQqOOOOO- 00.' 


GAIN ADOuSTHENT is HA'LVEOVS-,''i''-Sr r 


SYSJEH^IGENVALUES'V.'.'. V -'-:' ‘ 

■ ‘~QL4S07203b-'Ot::'~--C.O'OflOC‘aOD^'flO-:,' ■ - .-I.-. ,• 







FIGUBE 111.5-4 (Cont'd) 
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INVERSE GRADIENT HATRIX 
1.703A993 



FIGUBE 111.5-4 (Coni' d) 
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6a.254^iJI .es =5., OCQjlvq.OIL.gQ 
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FIGURE III. 5-4 (Cont'd) 
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FIGURE III. 5-^ (Cont»d) 



Average 

Cost 


2.0 


1.6 


1.2 


.8 


,h 


0 


FIGURE III. 5 





77 


on the computations. The systoa equations -were 

x^(t) 0 1 0 ' ’ x^(t)' ’ 0 ' 

X2(t) = 0 0 1 Xg(t) + 0 u(t) (III. 5 - 7 ) 

x,(t) 670 x_(t) 1 

■with the cost ftmctional 

[■30 o1 1 fioo" 

j == x^(i) 0 2 0 x(i) -i- J x~(r) 0 2 0 x(r) + 0.1 u^(r) ar 

_o 0 ij 0 ® 3 

( 111 . 5 - 8 ) 

It -was assumed that only the first state -was available for feedback so that 

u(t) = - x^(t) (III. 5 - 9 ) 

The computer • print-out for a starting gain of aero is given in Figure 
111,5-6 and is con^jared -with an I.C. in Figure III, 5 - 7 . There is an 
offset in Figure III. 5-7 be'tween the average costs produced by ISAFT and 
those produced by the I.C, This is due to the difficulty in approximating 
a uniform distribution on the surface of a sphere hy a finite ntmiber of 
discrete points. (This difficulty increases rapidly with the order of 
the system). It is nevertheless clear that ISAFT achieved a Tm'-nimuni 
average cost. 

Ill . 6 Conaments on the Finite Time Problem 

Ihe objective of the I.S.A, approach is to determine the feedback 
gains gi-ving the minimum value of the gain function GF> which is the 
expected value of the cost functional. The I.S.A. algorithm presented in 
Section 111.3^ however, is designed only to determine feedback gains 
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satisfying certain necessary conditions. The necessary conditions are 
that the partial derivatives of the gain function -with respect to each 
feedback gain should be zero. These necessary conditions may be satisfied 
by local maxima, minima or points of inflection. With a multidimensional 
problem points of inflection seldom occur in practice. The convergence 
factor u built into the digital computer program ISAPT prevents con- 
vergence to a local maximum. There still remains the problem of convergence 
to a local rather than a global minimum. Wo solution to this problem was 
foxind. It is felt, however, that it should be feasible to guarantee the 
existence of a single (global) minimum for suitable system and cost matrices, 
and useful work coirld be done in this area. 

The digital computer program ISAFT computes the various required 
functions of the state transition matrix (see (ill. 3-33)) through (ill. 3-38)) 
analytically. It was felt that this best illustrated the procedure. For 
high order systems however, it should be more economical to evaluate the 
state transition matrix at a number of discrete time instants and then 
evaluate the required functions by numerical integration. 
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CHAPTER IV 

SUMMARY AltD CONCLUSIONS 

This report has considered the determination of optimal time-invariant 
output-feedback controllers for linear dynamic systems with quadratic cost 
functionals. The need for such controllers in practical engineering was 
illustrated, and the contributions of earlier researchers in the field 
were reviewed. 

It was determined that deficiencies existed in two areas. For the 

case where optimization was to take place over the semi-infinite time 
19 20 

interval Levine ^ and Cassidy had each derived a suitable computational 

algorithm. Both algorithms, however, required initialization by suitably 

stabilizing feedback gains and neither author gave a method for determination 

of such gains. For the case of optimization over a finite time interval, 

no satisfactory existing techniques were uncovered. 

The problem of the determination of stabilizing feedback gains was 

approached via a gradient technique. The technique evolved from an 

27 

eigenvalue sensitivity relationship given in Fadeev and Fadeeva . It was 

mechanized in the digital conputer program GRADGN (Appendix 1) and provides 

a practical method for determination of local stability maxima in feedback 

gain space. It is felt that this gradient technique should complement the 

existing optimization algorithms for the semi-infinite time interval problem. 

A new technique was derived for the finite time Interval problem. The 

technique is based on the Initial State Averaging concept, previously used 

19 

for somewhat different problems by Levine and by Klelnman, Fortmann, and 

pg 

Athens . A computational algorithm was derived and is incorporated in the 
digital computer program ISAFT (Appendix II). The algorithm satisfies a 
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set of necessary conditions "by Rewton-Eaphson iteration^ using their 
gradient with respect to the feedhack gains. This is equivalent to 
minizization of the escpected value of the cost functional "by the laethod 
of second variations. 

The contributions made by this report should aid in the search for 
practical optimal controllers. 

Several outstanding problems remain. T'oremost among these is the 
question of the sufficiency of the solutions obtained by the optimization 
technique described above. The .technique was designed to determine local 
minima of the expected value of the cost functional. An examination of 
the convexity of this quantity in feedback gain space rai^t uncover 
conditions ensuring a single (global) minimum. The gain initialization 
technique is similarly local and would also benefit by extension to a 
global technique. 

The computational methods used in the digital eonputer program ISAPT 
were designed to illustrate the theory. They are satisfactory for low 
order systems, but not for high order systems. It is felt that relatively 
sinple modifications to ISAPT should eliminate this deficiency. 
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APPENDIX I 


This Appendix co^ri’ses 9, listing of the gain initialization digital 
computer program, GRADGN. 



1 


/JUB 


9H 


‘•257 KCf'kt NiN»LlMES = 5;> 

DI -icKSlDM H ( J ) »XK (3 > ) » X I ( 3 I , VK ( 3..i) , V I ( ?'>) , lAOKO 1,2 1, 

1. _ .. ._liS4ufJXaiU3i:L«JiSAI!li33,3:jJ,GR.AC(3j-.3JJ,.VB.'J13lU,JiLia3.Jl^ 

2 Gl3v,30) 

2 COMMON AHATOv, 30 ) , CHAT 1 30,30 ) , AHAT t 30, 30 > , AAAA 1 90 J I , ASOR 1 30, 30 > , 

L RRt3m_,KJI3'JJ.IAliAL3DJ_ .. ... _ _ >. 

1 kPAL K(3.',3J) 

- !..> FORMAT (71 IJ) 

. . -2^ , „-£JHajJ-JAEI3..4J 

o 25 FnRHAT (2<>AA1 

7 30 FORMAT (7F10,A-) 

r, 35 .format {’1',TA.' STATES' . AX . « CONTROI S ■ . AX > ' FgFOR ACK ?; « . AX . ■ fS AINS'.,. 

1 ‘X, ' i;j=LR' ,AX, • ISTOP' ) 

) 36 FO-IMAT (/7T3, 'STABILITV INCRfcASE INCREMENT =',F20.8) 

I - 37, ...FORMAT .lZZIl. Li _ M_L NI MiLM ST AB [ L l.Ty_..R_£QUJ&.F[J .FOR . .. 

1) 

II 4u FORMAT (//T3, 'STSTEM MATRIX AMAT*) 

12 F0RMAT{//J3^.«C0NT R0L MATRI X BHAT' 1 .. . . . . 

;3 jtj FORMAT {//T2, 'MAXIMUM REAL PART UF ROOTS'! 

lA 6.- format (//T3, «E JGVEC ERROR MESSAGES') 

I .. . FORMAT ( T3.’S„1 = » ►E10.4,IOX . ■ 1 T ER=. ' . 1 5 . IPX . ' C 1 F= ' . £ 1 0 . A 1 .... 

lo 7v’ FORMAT ( //T3 , ' 6 IGBNVECTOftS CORRES TO MRP EIGENVALUE') 

17 75 FORMAT (T6,'R0W REAL PART' ,5X, ' ROW IHA5 PART* , 5X, ' COL REAL PART', 

. _ l_5Xa'CDL..lHAG.J‘.ABJJ.l . 

II a.J FORMAT (//T 5, 'GRADIENT MATRIX') 

n 82 FORMAT (' 1 ' ,Tj, 'PROGRAM TERHINATIO-J DUE TO EXCESSIVE NUHpER OF ITE 

IKAT irvs') 

2. 8A FORMAT ( ' 1' ,T3, ' SPEC IF lED STABILITY HAS BEEN ATTAINED'! 

21 86 FORMAT ( ' 1 ' ,T3, "GAIN SEARCH PROCEEDS ALONG A NEW GRADIENT. ITERATI 


2.1 

2 J 
2 A 
25 

. 26 

2 7 


',. 301030 . 

87 FORMAT (//T3, 'GAINS TOO FAR FROM STARTING PtJiNT. COMPUTE NcW GRADI 
lENT') 

88 FORMA T (//T3., 'STABILITY INCREASE STEP SI Z E . F 2 J . 8 ) 

VO FORMAT (//T3,'STEP SIZE IS DOUBLED* 1 

92 FORMAT (//T3, 'CONVEX CURVE FIT TO STABILITY LOCUS'), 


95 FORMAT (//T3, 'CURVE FIT NOT SUCCESSFUL. RETU.RN TO PREVIOUS BEST GA 
IINS') 


27 

3 ' 
3i 


95 FORMAT ( 7/ T3_. 'STEP S I Z.E_T.CS_ L ARGE. COHP.UJ E NEW GRADIENt_A.t_PREVIO U 

IS BEST GAINS' ) 

96 FORMAT' t • 1 ' , T3 , ' STAB ILI TY IMPROVEMENT RATE TOO SLOW, PROGRAM TERM! 

l.NAT £.0,..'J : , 

97 FORMAT (//T3,'LAST RESULTS SHOW HOST STABLE CONDITION FOUND'.) 

98 FORMAT (//TB.'STEP SIZE TOO LARGE. PROGRAM CONTINUES WITH REDUCED 


ISTeP SIZE') — 

C 

C NS NUMBER' OF SYSTEM STATES 

„ Q NC V NUMBER "OF CONRTROLS 

C NF = NUMBER OF FEEDBACK STATES 

C IGAINS = ') - INITIAL FEEDBACK GAINS ARE ZERO 

C = 1 - INITIAL FEgOBAOK GAINS AR E READ (N 

C IDELR = 0 - STABILITY INCREMENT = U.'l 

C = I - STABILITY ENCREHENT IS READ IN AS OELR 


C ISTOP = 

C 



1 - PROGRAM TERMINATES AT STABILITY SPECIFIED BY STOPS 





32 

33 


C 


READ (1, 13) lOELK, ISTCIP 

HBIXS. 

34 WRITE (3,IJ) MS,NC,NF,IGAINS,IDELR,IST()P 

C 

C nSLR.CONTR GLS STABILITY INCRE ASE ^ICREHEHT. J.. . . 

33 IF IIOELR) 12-J,123tU:i 

36 llo READ II. 3H) O.M.R : 

37 GO TO 130 

38 120 DELR=0.1 

39 130 DELR1=0ELR 

4-1 WRITE (3,36) uELR 

C 

C STOPR CONTROLS STAtilLlTY RFOUIRFD F UR PROG RAM TFRH I NAT IU3 

C 

41 IF nSTOP) 150,150,140 

42 _ 140 READ 11,30 1 STOPR 1 

43 GO rn i6j 

44 150 STJPR=-1.9S6 

45 loo WRITE (3,37) jTOPR 

46- DO 170 1=1, NS ' 

47 . 00 170 J=1,NC 

48 170 Kn,J)=0, 

c 

L SYSTEM DATA MATRICES READ IN RY ROWS 

C AHAT IS SYSTlM MATRIX 

C BHAT IS CONTROL MATRIX 

C 

^ REAP (1,301 (1AHATII.J),J=1.NS1. I=;1 iMS.1 , 

5j WRITE 13,40) 

51 WRITE 13, 2J) ( I AMATI I ,J) , J=1,NS ), r=l,NS) 

52 READ (1,30) t ( BMAT <1 , J 1 , J= 1 ,NC 1 . 1 = 1, NS'i 

53 WRITE (3,45) 

54 WRITE (3,20) ( ( BMAT ( I , J ) „J=1 , NC I , 1= 1, NS ) 

c : 

C INITIALIZATION OF GAIN MATRIX K, 

C 

55 IF (IGATNSl 1 9 J . 1 9.3 ■ 1 8 3 

56 180 READ (1,30J ( IKI I , J » , J=1 ,NC ) , I=1,NF ) 

57 190 KDUNT=1 

58 • ROOTT=0, ; 

C 

C SUBROUTINE STAB COMPUTES THE SYSTEM STABILITY 

C 

59 195 CALL STAB ( NS ,NC ,NF,K) ' 

60 198 IFIRRUI-STOPRJ 999,200,200 

61 200 KaUNT=KCUNT-H 

62 IF (KOUNT-U) 2iO,2Ui,ll..<j 
C 

C GAIN SEARCH PROCEEDS ALONG A NEW GR AD I EN T 

C 

63 210 RO0TR=RR(i) 

64 ftOOTI=Rl(l 1 • __ 

WRITE (3,53) 


65 



60 

67 

63 


WRITE (3,2J) <snGTR 
RT6ST=ROOTT 

ROOTT^^ROOTR 

69 K1=K0UNT-1 

70 WRITE (3,86) Ki 

£. 

C CQMPUT4TI0N OF EIGENVECTORS 

C 

7i CALL MSQ(NS) 


72 CALL tIGVEC t3,AHAT,ASQR^W,IR0W,XR,XI ,VR,VI,ROQTR,RaOTI,NS,3D,0, 

1SW1,ITER,DIF,2) 

C 



C 

C 

C 

SWl = J for' EXACT EIGENVALUE AND NO RCUNO-OFF ERROR 

ITER = NUMBER OF ITERATIONS USED TO pjND EIGENVECTORS. 

I TER = 15. 


c 

Q 

DIF = LARGEST CHANGE IN ANY EIGENVECTOR COMPONENT AT FJNAL ITER. 

. JJL 


WRITE (3.60) 



74 


WRITE (3,65) SWl , ITER, niF ' 



73 


WRITE (3,7 •() 


'' 

7o 


WRITE (3.75) 



77 

c 

C 

WRITE (3,20) (VR(l:),VI(I),XR(I),XIin,I»l,NS) 
NORMALISE FIGFNVFCTORS INNFR PRODUCT 




c 




•7<i 


V£CMGR=0. 



79 


VECM.GI = 3. 



80 


DO 2A0 1=1, NS 



81 


VECMGR=VECMGR + VR 1 1 ) *XR( I )~VI ( I ) *X 1 1 1 ) 

\ 


„ „ 82_ _ 

2<V0 

VECMGI=VECHGI+VR[ I ) *X I ( I ) *91 ( 1 1 *XR ( T ) 



8 3 


VEC«GS=VECMGR»VECMGR+yECMG[*VECMGI 



89 


DO 253 1=1, NS 



65 


VRN(.I) = (VR( I) «V£CMGR4VI ( I 1»VECMGI)/VECMGS _ 



86 

250 

C 

VIN(I) = (Vm) *VECHGR-VR 1 1 ) *VECMGI ) / VECH6S 




£ COMPUTE KR API ENT MATRIX 

C 


87 


OG 3C0 J = 1,'1F 


88 


on 300 L=1,NC 


89 


GRADR( J,L)=0. 


90 


GRADnJ,L)=0. 


92 


on ?90 T=1 .NS 


97 


GRADR( J,L)=GRADR(J,L)-VRN{n»BMAT( I, LI 


93 

290 

GRADI ( J,L)=GRAOI{J,Ll+VINn )«BMAT(I,L1 


94 

300 

GRAD ( J.L) !?GRAuR(.I.L) «XR ( J ) VgRAO I ( 1. 1 ) ♦Xf^.l ) 


95 


WRITE (3,80) 


96 


WRITE (3,20) {(6RA0(I,l}rJ=l«NC), I=I,NF] 


97 


WRITE (3.88) DELR 


9b 


GR0S0=0. 


99 


DO 320 1=1, NF 


l-n-j 


DO 320 J=1.NC 


lOi 

320 

GR0SQ=GRD$0*GRAO(I,J)*GRAD(I,J) ' ' - - 


102 


DELKI=-DELR/GR0SQ - 


1-J3 


DO 340 1=1. NF 


104 


DO 34C J=i,NC 

' 



i;t5 

K>6 

1'j7 

34U 

vJ (I ,J)=K(T , J) 

KU , J)=Gt I , J) fOeiHluGslArH I , J) 
pnijri=RooTK 

108 

109 

110 

390 

CALL SrA6tNS,iVC,NF,K) 

If tRRi D-STOPRI 999,390^390 
RCUT2=RRf 1 » - ... 

111 

112 

113 


KnU'lTl = l 
K.0U‘JT2 = 1 

IF iRonrz-Roam v.i i.4n.).6-ij 

114 

400 

IFnROOTl-ROOT2)/DELR-,95) 420,420,500 

IIS 

420 

DELK2=,5*DELR*DELKl/IROOT2-ROOTl-FDcLR) 

116 


WRITE (3.92) . 

117 


DO 44C I=1,NF 

118 


DO 440 J=1,NC 

119 

440 

KfT*J)=G( K?*fiRAn( T«.i) 

120 


CALL STAa(NS,NC,NF,K) 

121 


IFlRRtU-STOPR) 999,999,460 

122 

460 

R00T3=RR(1) • 

123 


lFtR00T3-RC0Ttl 470,47), 480 

124 

470 

IF(K0U'JT-2) 4?5,475,472 

125 

472 

IF(RTEST-RCDT3r-DFLRl 1 12Du.475.475 

' 126 

475 

DELR=.5*(ROOTi-RCOT3) 

. 127 


GO TO 20'J 

128 • 

480 

IFIKOUMT-2) 485.485.482 

129 

482 

lF(RTEST-RC0T2-DELRi) 483,485,485 

13J 

483 

DO 434 1=1 ,NF 

131 


DO 484 J=1 ,NC 

132 

484 

,J)^DELK1»6RAD(I,J> 

133 


WRITE (3,93) 

134. 


CALL. STA8(NS.NC,NF.K) 

135 


GO TO 1200 

136 

485 

DELR = .5*(R0(lTi-RC0T21 

137 


rtRlT.E_.(3,94) 

138 


DO 490 1=1 ,NF ■ ■ . 

139 - 


DO 490 J=1,NC 

140 

490 

Ktt,J)=G( I .J)+DELK1*GRAD( I. J1 

141 


GO TO 195 

142 

500 

0ELK2=2.*0cLKi 

143 


KnUNTl=KnU-NTI + l 

r-44 


IF(KOUNTl-ll) 505,505,502 

145 

502 

DELR=.5*(RQOTi-RCOT2) 

_146 


WRITE (3,87) 

147 


GO TO 200 

148 

505 

WRITE (3,93) 

149 


OD 510 t=l.NF 

ISO ■ 

, 

DO 510 J®1,NC 

i5i 

510 

Kn,J)=Gn,J»*DELJ(2*GRA0(I,J} 

152 

* 

CALL STAB(NS.NC.WF,K) 

153 


IF{RR(1)-ST0Pk) 999,999,520 

154 

520 

RnaT3=RR(l) 

155 


IF (RUUT3-ROOI2) 530.531.58j 

, 156 

530 

ROOT2=RQOT3 

' 157, 


0ELK1=0ELK2 

-'■158 . 

- « 

60 TO 400 

159 

580 

IF(KCU(^T-2) 585,585,582 



96 


16j 582 IF(RTEST-RR1JT2-DELR1 J 583,585,585 

161 583 on 58^, 1=1, NF 

Lii2 PQ 5 fl <t...J = l ,r< C 

163 58<V Kn,J)=G(I,J)-fOELKl»GRAD(t,J) 

164 WRITE (3,93) 

1&5 CALL STAB(NS.NC.NF.K) 

166 GO TO 12UO 

167 585. OELR=,5»(ROnTl-RCOT2l 

LM WRITE (3.951 

169 DO 590 1=1, NF 

170 DO S90 J=1,NC 

.Ui 53J} K(T ,J1=GI I,J)*OELKl*GRAD(t.J) 

172 GO TO 195 

173 630 KOUi'JT2 = KOUNT2*-l 

12A IF(.C0im2-61 0 35.605. im 

175 605 OeLK2=.5*DELR»OELKl/(ROOT2-RQOTl+OELR) 

176 WRITE 13,92) 

111 00 620_ I = 1,NF 

178 DO 620 J=1,NC 

171 62J K II , J)=G( I , J)+DELK2*GRAD(I , J) 

. -iPv cmj._jiiAPj N s , n c , N P. ,8; ) , 

181 IFIRR(l)-STOPR) 999,999,640 

182 ■ 640 R00T3=RR(1) 

-L8J tF(R00T3-RCDTl) 680.680.6^0 

184 660 0ELK1=DELK2 

185 WRIT? 13,98) 

1.86 RnnT?=RnnT3 

187 GO TO 600 

188 680 IE(K0UNT-21 635,685,682 ' 

189 682 IF(RTEST-R00T3-D£LR1) 1200,685.685 

19,) 685,. DELR = .5»IROOTi-RCOT3) 

191 ' GO TO ZZO 

192 999 WKITF (3.84) 

193 GO TO 2000 

194 1100 WRITE (3,82) 

IS5 GO TO .2000 ^ 

196 120-' WRITE (3,96) 

197 WRITE- (3,97) 

198 20Qj CONTINUE 

199 STOP 

20a END 


I 





2'Jl 

232 


SUBROUTINe ST^6(NS iNCtNF.Kl 

COMMON AMAT(i^,30) ,c;MAT(3 i , 3 J 1 t AHAT ( ^.i , 3 1 ) , A A AA ( 9 J J ) , ASU-<( 3,»,33>t 

203 REAi~Kt30i30> 

20A 10 FORMAT (//T3,*GAIN MATRIX K' ) 

205. 20 FORMAT (AE18.7I 

2j6 30 format (//T2, 'ROOTS' ,5X, 'RIAL P ART ■ , 1 IX i ' IMAG. PART*) 

237 4i/ format (2E2),o) 

203 ■ 1A = HS _ 

209 WRITE 

213 WRITE I3,2i| I I K I I , J I » J=.l ,NC ) , 1 = 1, NF 1 

211 CALL AHHTINS. NC.K) , 

212 CALL VECT(NS) 

213 CALL HSRGlNS.AAAAflA) 

21A C ALL A T EtGINS.AA AA .RR.R I.I AN A. lAI 

215 WRITE (3,30) 

216 WRITE (3,A-.)I (RR(I),Rim,,I = l,NS) 

217 CALL HAXR TINS) . 

RETURN 
ENO 


2ia 

219 



22 


subroutine AMiiT (ns tNCi K) 


C 

C rOMPUTrS AHAT AMAT - BHAT.K ( T RAMSPUSP I 

C 

221 COMMON AHAT(3o,30),BHAT(3i0,3O),AHAT(30,30) , AAAA(900»,ASQR(30,30) . 

I RR n q.) .^ulao)«IA f ^A o o) 

222 REAL K( 3^1,3 )l 

223 DO ICO 1 = 1, 'IS 

22 ± on ICJ J=I.N S : 

225 AHATd ,J)=iAMATn , J) 

226 00 100 L-1,NC 

227 100 AHAT(l.JlnAH A T(I.Jl-BMAT(l.Ll»KtJ ^Li _ 

22fc RETURN 

229 END 



99 





238 


239 

240 

241 
-i*2 

243 

244 
.245 
246 


subroutine HSQ(NS) 

C 

C . COMPUTES JVSOR AHAT*AHAT 
C 

COHflON AHAT00.30 I ,6HAT (30*30), AHAT( 30,30 )• AAAA(900) ,ASQR( 3O,>30> , 
1 Rft(30),R[<30l*IANA{30) 

DO TOO I-1,NS 
00 100 J*1,NS 
..AS(»iJl*Jl30> .. 

00 100 K»1,NS 

100 ASQRM* J) 3 ASQRd ,J l-l-AHATl I ,K ) >AHAT( K, J ) 

. . RETURN . . 

END 



249 SUBROUTINe HSQG(N»A,IA) 

C 

C- CnN V£R_T.S A T0...upPgR HES S gNBrWG . FORM 

c 

250 OOUBLi: PRECISION DABS t OFLOAT «DS IGNi OBLE i DEXP, 0LO6,0L0&l0« OATAN 

l.DSIN.D[:OS.DSgRT.>DTANH.QHOO.DHAXlwDHINl 

251 DOUBLE PRECISION S 

252 DIMENSION A(9.‘JI 

£12 Uzii 

254 NIA=L*IA 

255 UA=NIA-IA 

216 20 IF(L-3) 360.40.40 

257 40 LIA=LIA-IA 

255 L1=L-1 

259 L2=L1-1 

260 • ISU0=LIA+L 

261 IPI\/=ISUB-IA 

262 PIV=ABS(AUPIvn 

263 IF(L-3) 90f99,5J 

264 53 M=IPIV-IA 

263 on 89 I^L.H.IA . 

266 T=ABS(Ann . • ■ 

267 IF(T-PIV)' 80,80,60 

268 60 IPn/=«I _ . 

269 P1V=T 

27J 8J continue 

271 90 IFIPtV) 10.J.323.1CI 

272 100 IF(PIV-ABS(A(ISUB)n 180,180,120 

273 120 H=IPIV-L 



285 - 160 A( J)=T 

2 86 1 8 j_D0 2iS-J.= k, UiLiJL&. 

287 2UJ A( i) =Atn /4(I jUBl 
28,i J=-IA 

289 DO 24J L=_1,L2 _ 

29i) J=J+IA 


291 . LJ=L+J 

292 DO 2 20 K=1 ,L1 

293 KJ=K+J 

294 KL=K+LIA 

293 2 2J A ( KJ ) =AJ_K_J )-A t LJl *A 1 KLl 

296 240 CONTINUE 

297 K=-IA 

298 DO 3.)0 

299 K=K+fA 
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314 


suBRQurrNe atcigi m,a,rr,ri,ijina, ia> 


C - COHPUT&S'ROOTS OF UPPER HFSSENBgRK MATRTX A 



DIMENSION A (900) ,RR (30 ) »RI (30 ) , PRR( 2) , PR I ( 3 ) , IANA(30) 
INTEGER P,P1tO 
B7=1.0E-8 



HAXIT=30 
N=M 
20 NL=N-1 


30 NP=N+1 
IT=0- 

00 40 T=1 



N1N=IN+N1 
N1N1=IM1+N1 
60. T=A(N1N1)-A(NN1 



65 T=U+V 

IF(ABS(T)-AMAXl(U».ABS(Vn»E6) 67,67,68 
67 T=0.0 



70 IF(U) 80,75,75 
75 RR(Nl)=t+V 
RRINt=II-V 



-358; 

' 3’59; 


GO TO 130 
100 IF(TI 120,110,110 
110 RR(NU=A(N1N1I 






i6‘^ 13.1 RI('n= I. .1 

i6o = 

3.6 7 r.n in i6:t 

36R UO RRlNl)=t 

369 RRtNJsU 

37i( M(NU°V ; 

371 RHN)=-V 

37.; 16J IF(N211280,12JJt18'» 

- 37 > l80.._VLN2.?Nm-I.A 

374 RHaD=RRlNl)*RRtNl)+RHNll*RI(NU 

375 ■ EPS=ElO*SORTtKMOD) 

3.7.6 IF( ABS(A ( N1N2I )-EPSI128a>12Bl).24u 

377 24o IF(ABS(A(fm) 1-E13*ABS(A(NM n 13:i J , 13 , 25J 

37i 25vl IF(ABS(PAN^-A(^llN2l )-ABStA(NlN2n*e6) 1240,124 1,2 61 

. 370 26J IFlAB.St PAN-A<i'i,Nlt )-ABS(AtNNl 1 > »F6 1 1 24Q . 1 24.0. 30.1 

38 ( 300 IFilT-MAXIT) 323,1240,1240 

381 320 J==l 

.382 P Q ..36D_JL=JU2 

383 K=MP-1 

38 A IFtA3StRR(K)-PRR(! ) 1 + ARS ( RI ( K ) -PR I I I) ) -D ELT A* ( ABS ( RR ( (O ) 

1 4ABSCRI tKl 1 1 1 3.4J.362,J6J 

385 340 J»J+[ 

386 360 CONTINUE 

3 37 r.n TO (443. 46J. 460.4BQ1.J 

388 44V R = 3..1 

389 s=3,n 

. .ATj aP— LQ...?J3 

391 460 J=N+2-J 

392 R=RR(JJ*RR( J) 

. 3J93 S=RR< JH-RR( J1 

39a go to 519 

395 483 R=RR(N) »RR (N1 )-Rl (M •RI (N1 ) 

396 S=RRINl+RR(N ll 

397 500 PAN*A{NN1) / 

398 ~^PAN1=AJN1N2) 

_3.95 &0-520 1 = 1. 2-^^ 

4Ji Kl=.Nt>-I . 

4'>1 PRRm=RR(Kl 

4j; ._..52 j PRl ( I)=RHK) 

405 P=NZ 

404 IF<N-3)630,60'J,525 

5?6 IP1=N1N? 


406 on 580 J=2,N2 

4,17 IPI = IPI-I4-1 

4 08 !F<ABS(A(IPin-EPS) 6Jj.601.530 

409 S30 IPIP»IPI+1A 

410 IPIP2=IPIP+IA 

._4JU n*A(IPTP>>[AI IPIP)-SI»AHPIP2l»AnPIP-Hl4R 

412 IF(D)54J,56J,j4U 

413 540 IF{ABSUnPI)«A(IPI1> + l) )*(ABSIA( IPIP)+A( I P I P 2+1 )-S ) +ABS 1 A t I P I P2 + 2 I 

1 1) -ABStPl»£PSJ_.62J.,.6.2iJ>560 

414 560 P*N1-J 

415 580 CONTINUE 

. 4 1b 600 Q=P 

417 GO TO 68j 



A18 

419 


620 Pl=P-l 
0 = P1 


'650 OO 66J 1=2, PI 
IPI=IPI-IA-1 

IF(ABS(A(IPH )-EPS)6B0.6B0, 

6oO 0=Q— 1 

680 I I = (P-1)»IA+P 

t)Q 122-> t=P.Ni ^ 

Ill=II-IA 

iiP=n+iA 

IF (I-P 1720,700,720 

70 ) IPI = II + 1 
!PIP=IIP+1 

Gl=A(n )»tA [I 1 l-<;i -fA{IIP)»J 
G2=A( IPI)»(A( 1PIP)+AII I )-S 
G3=A(IPn»AnPIP+l) 

A(lPt*ll=O. 0 

GO TO 78j 
72.J G1=A(1U) 

G2=A(IIItI) 

IF! I-N2I74 J,74'J,76l 
740 G3=A(1I1+2J 

GO TO 780 

76 j G3 = J • 0 

78 J CAP = SutlT(Gl*r.i+G2»G2+G-3»C;3 

IF(CAP>a.-n t86 ->,60 0 

800 IF(G1)820,840,840 
820 CAP=-CAP 

640 T=G1*CA P 

PS 11=1,2/ T 
PSI2=G3/T 

ALPHA=2. .1/(1. itPSIl»PSlMP 
CO TO 880 
860 ALPHA=2.0 


PSi2=-j.O 

8bJ IF(I-O) 9-j 1,96 

9tlJ IF t I-Pl 92,1.94 .,92 i 

920 Anil)=-CAP 
GO TO 960 

940 At rtn=-Atiti 1 

96j IJ=II 

00 114'J J=1,N 

T=PSU*Anj-H) 

IFn-Nl)980,lw00,lQU9 
980 IP2J=IJ+2 

• t=T-«-PSI2»A.MP2J ) 

1000 eTA=ALPHA*(T+Anj) ) 

A I I J)=A( IJ l-= lA . 

AHJ4l)=A[IJ4-il-PSIl»PTA 

IF (I-Nl 11920, i;i40,lu40 
1020 A(IP2J>=A(IP2J1-PS12*ETA 

1040 IJ=IJ*IA 

IF( I-,'J1)1J39,1j6C, 1.160 




47 3 

lo60 K=N 


474 

GO Tl) 1190 


475 

ItlR.) K = t+P 


476 

IlOO' IP=HP-I 


477 

00 U80 J=50fK 


478 

JIP=1P+J 


476 

JI=JIP-IA 


48. 

T=PS11*A( JIP) 


481 

IFn-Nl)112'J,I140,.114:i 


482 

1120 JIP2=JIP+IA 


483 

T:*T4PSI2*A(JIP21 


484 

1140 £TA=Al PHA* { T+A( JTJ 1 


485 

A(JI 1=A( JM-C(A 


48o 

At JIP) = At jip)--:ta»psii 


487 

IF( I-N11116:i,xlflO*llP'j 


483 

1160 At JtP2l=A( JIP21-ETA*PSI2 


489 

1160 CONTINUE 


49 "J 

tPt I-N2)1200, 1220*1220 


491 

1200 JI=II+3 


492 

J!P=JI+IA 


493 

JtP2-JIP+IA 


494 

ETA=ALPHA*PSI2*At JIP2) 


495 

At JI 1=-ETA 


496 

AlJIPI=-FTA«PSn 


497 

At JiP21 =At JIP^)-£TA»PSI2 


498 

1220 I1=I1P+1 


499 

IT=IT+1 


50J 

GO TO 60 


501 

1240 IFtABStAtNNll)-ABSlAtNlN2m 1300,1280,1280 


•5i)2 

1280 !ANAlN)=n 

'* 

5v'3 

IANAtNl)=2 

' 

504 

N=N2 


a'lD 

IF(N2)'l4-jJ.14yJ,20 


506 

1300 RRtN)=AtNN) 


50 7 

RT(N)*?0.0 


508 

IANA1N)=1 


599 

■ IFtNl) 140'J,14jl),132a 


Sly 

1320 N=N1 


jl 1 

GO TO 20 


512 

1400 RETURN 


513 

END 




5i4 


SUBROUriME MAXI^TCNS) 


C 

£_ COMPUTES _RQOT HAVING MAXIMUM i< £.^L PART 

C 

515 common AMAT(30,30) tBKAT(3!),30)f AHATI30,301,AAAA(900lfASOR(30t.30) 

1 Rdt.am iRI t30Li.UNA.t.iQJ 

516 DO IJO 1=2, NS 

517 IF UR( IJ-RRI in -Su,! j1,l jj 

518 5u RIlll=R[m 

519 RRm=RR(I) 

520 100 CONTINUE 

521 RETURN 

522 END 



lo8 


523 SUBROUTINE EIGVECUVC, *, 6» XR, XI, VR, VI, ROOTRE,- ESVI. 

1 ‘r'OOTIE, NE, NMAX, T2, SHI, COUNTS', ERR.MMHI . ; . SSVl, 

C SUBROUTI^NS TO EINP' THE EI^EN VECTORS , QF A NON-SYMU^RiC HATftly SAXJL 




c -- 


2 FINO 0/{tY, THE TRANSPOSER EIGENVECTORS (AT V « LAH80A VIESVl 

3 FIND 'MTH .TYPES OF EIGE&VECTORS. ‘ ESVli 

dimension A(3Q,3P1 rB?3D,3Q» ,W(30rAi>XR(3pi .XI f3ni.VR( 3Ql.tf 1130.1^ 



563 


622 1F(IVC2) 623, 625, 623 




esvi 
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570 626 VULC JsO.O' ' 

571 GO TO '5li 
JL 






- ; I2';f=t IROHfIfZl 

' ■ \ 1 > • H < I . U »RGOT 1 

- ortVfil i' .1 = 1 . N 



= »( I »3 MRQ.OTI 
• DO 61’3 J 3 ItN 
613 -VTTi) = VnU 



616 DO 618 I =1 1. N 
•. ■;?«(!) = -W(I,2) 
-PD '617 J .=i I, N 



597 

619 DO 621 1 =? 1, N 

598 

-.'vIki) =-Hn,Aj 

599 

DO 62&-J.3 1. N 



SEARCH VECTORS FDR LARGEST ELEHENT ANC NORMAL U6, 



606 

607 

606 


628 AMAX = TEHP 
12 = L 

629 CONTINUE 


609.'--/v;;:^a-.(?r:= VR'(T2 ) /AHAX 
- -il {'I2 ) /AMAX 
63g-C =; 1. N- 


612 


TEMP = VIIL) 


ESV 

£sy 

ESY 

esY 

£SY 

EGY 

C-SY 

ESY 

cSY 

ESY 

CSY 

cSY 

cSY 

ESY 

CSY 

ESY 

ESY 

ESV 

ESY 

ESY 

ESY 

ESY 

ESY 

ESY 

ESY 

ESY 

ESY 

ESY 

.ESY 

ESY 

ESY 

ESY 

ESY 

cSY 

-;sY 

ESY 

ESY 

ESY 

ESY 
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yi‘(U = VRtL)»C2 ♦ TENP*C1 
630 VR?,L> = VR(L»»C1 - T6MP»C2 



636 XR,a) = XR(U»C1 - TEMP»C2 
IFtCOUNT’-.eQ. 1) GO TO ‘6<»6 
DO 637 LL = 1, N 



638 IFJCOUNT ,E0', IJ GO'TO 646 --‘i:''i'.. ' 
IFtCFRR- -CP* CH’ Tn 







111 - 



ESV3 

esy; 



•-.668::V;.7/'652-.IFfrrCCK 680,-. 685-,_r680 ' r: 



,:67ls>;.K.682:xrt,U' , ,.. , - 

f.67g .; -I "FlI VC 1 ). 6 83 , 5 lV^*:683sH4.^ ••; \ 

ri673~:-''.4-683 OQ. LVa ' iC 'M ~ -"->-'W., ^ 



,,6'8l«-,!;i4657r^.i ( C V : V; *' ' ^ : 


kESY, 
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708 

i:70g 
■ rici 


CERft = AHA'Xl{CERR> ABS.(V,ldL)' M,( LU.l ) )'} , - 
wautiJ^vuLcr • " 

665 VRn.L) = H(LL,l) • - . i . . 

* ' * 

.. 'v.-Es^ 

S 71-1- ' 
f 7.12 
1 7,13 

-- 

GO TO 638 ■ , , 

66.8 EPUVC2) -669,- ’6:7l> 669 ’ ' VJ - r 

669 DO 670 L- 1 .. /J - - - • -.-V ■ ~ 


. ,.,.'6SY 

ESY 

-i 

i 

fTnT 

} 7.15- 


67j3 XHU = 0.0 . . . 

IFUVCl) 671, 70,. -671 ' ■ 

671 DO 677 1- = T. M - ' - . ■ 

■ - 

, -V’ ’JESY 

ESY, 

• ‘ ‘ ’ ESlt 

? 717 
; 713 
i 7-W 


6'72 VKU) =0.0 

70 RETURM ' ' “ . ■ - 

673 IFnvC25 674. 5D7., 674 . * 

* , * 

- ... - ESYi‘ 

- - ■ '.ESY'- 

720 

721 
. 722 


.674 DO 675 I a, N 

12 = JROHn',2) ' ■■ ■■ • 

675 XriI2)' = XR(T'l - 

' > • 

, - ■ 'ESYl ■ 

, , .'ESY|; 

• : . _ _ . ESV. , 

} 723 

i 

0 

GO TO 500 

BACK SUBSTITUTION SECTION. • - ’• 


K- ESYi 

.esy;’ 

1 724 
' 725 

C 

499 IF(-EVC2) 500, 502^ 500 ' ’ 

500 DO 501 I s 7, ,'I . ‘ . 

■ . ■ 

ESY; 

. = , ESY'; 

• .FSV’ 

• -726 
1 727, 
728 


IX = I - 1 

- ■ >■ '• DD' 501- J =5 1,- M ' 

501 XHU =. xiMi - Rn..n*xi{.n - - ■ . - 

r ‘ 

.. -'.ESY 

■ - - "--Esv: 

"F!S¥V 

729 

.730 


511 IFUVCl) 502, 51-4. 502 ] ’ . . , 

502 DD 510 I = X, N ‘ - ’ - • ' 

n = 'l - 1 ■ - ' ' 


. ' - esy; 

ESYj. 

■ • ■ - E'<:v’ 

732 
733- ' 
734* 

- 

Tf(M) 503, 505,. 603 
603 00' 504-J’ = IV,. U' ' ■ 

504 vnn = VT.u-r- BKr.D.vifin': ... -- 

' • " f *' 

1 /V_ 

.-;,U --.‘v • ESYli 
‘ - ■.;• -ES7:l 

735 

736 

737 

1 

. -IFMCC) 505, 506, 505 . \ . 

505 IF(BUtJ)) 506, 507, 506' , ' . - 

506 Vim 4 Vim/BU.Tl . 

J ‘ 

■ ESY; 
ESV; 
FSV 

7'38 
» >39-.- 
1 .740 


■ ■ GO'XO 510 ■ • ■ ‘ . 

507 JFCVI-m) V508:, .509, 508 " '' ' ’ ■ " 

508’ vn-n ■= VMii*i>inp+i5 ■ -> . ■ 


,Vr-'ESV: 
• ! -' .■ : •- PSV.' 

741 


.■ GO TO 510 ' . ' 


esy: 

742 

743 


•509 Vim = 1.0 - ■ ' ■ - 

.5.10 CONTINUE- , 

' 

'ESY. 

FCV' 

j 744 
>■745= 
t -746 '■ 


;■ IPUVC2). 5i4i 525f; 514* ., . , 

-;514 DO 522 I-'=-'l-y--N. ‘ " - - ' ' 

IR - NEI 1 --I . - 


.- --.'esy; 

. - i.' , .... .esy: 
~ 'V • ' F.<:v: 

747 
-.74-8 
, 7_4a_. 


IFII - 1) 515, Sl'T, 515 . • - 

515 12 '= IR + 1 • ' - ■ ■ 

DO 516 . „U = 12. '-N ’ . . • . ^ - ■' 

• , , 

•- ESYJ 

'• • esy: 

■ PSV< 

Nso“ 

• 751. • 
? 752 


5X6. XIURl = XI UftJ BtJR,J).«XUa), . ' ■ . . 

IFUCcj- 517.,'5L8,“’517 • " ‘ ‘ ‘ 

517 IFIBMR.IRU’ 51».' 5l9’, -518 


-ESY-: 
,L esy: 

"psy-. 

f 753’ 

1 7.54 
f 755 ’ ' 


5l'8XUIR) = XIURJ/BURURJ ' ‘ - 

■GO TO 522 

519 IFIXt.UR).) 520. -521, 520 • . . T-'’ - 


"• • es¥l 

; ESYj 
• fsy' 

t 756• 

■-757.-r 

1758'. 


'520'XlflR)' - .XinR)*l-,0E+15 . . 

• , -GO TO 522’'' •' ■ ' ■ ' 

521. xmR»>= i.n ... 

'' ' *’*• ; 

.^^•,■::X .. .ESVf. 
.-:•“ .-n -Fsv< 

,-759 


522-;C0NTINUE ■ . . ' - . 

V* •' . , 

- - - : »ESy| 
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7j9 IFdCC .L5. lofti GC TiJ 71, 
1-.-5 I, 1 1 i51 

■ IvjSl . HRtTE(i n3. 1..I2. 1 ICC 

COUNTE = 0 
RETURN 


71 J on 711 LL - 1 , N 
12 = IROw (LL.l) 

711 IRUh'lI2.21 = LL 

IF(ROOTI) 6^17, 652, 607 
1052 FORMATI ///23H WARNING •*«»*» 


SUBROUTINE EIGVEC HAS 6SY 


1 F OUND AN E IG EnVAL UE OF , A£P A EE NT. JiULTllPJ.1 CUltl*. 

1 I4,/23X,* 

2Ge!'ivEcroft(si continues at user s option*//! 

iJl F0RHAT138H IMOk E THAN 13 LOOPS FUR ErGENVECTOR 

2 14H DIFFERENCE OF, Eli. 4) 

102 F0RHAT(16H0»*»»WARNING«**» , 14, 71H ZEROS ON 
I MATRIX. CHECK FOR MULTIPLE E IGENVALUES. /20X, 


ULTLlPJ^lCUlClj ESy 

I4,/23X,* COMPUTATION OF EIESy 

OPTION*//! 6SV 

ESY 

71H ZERCS ON DIAGONAL OF FACTOREDESY 

A1V,AL!JES./2QX, .ESY 


2* SUBROUTINE lIGVEC WILL NOT PERFORM COMPUTATION FOR THIS EIGENVECESY 
3TOR •//) ESY 

END 
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APPENDIX II 


This Appendix comprises a listing of the I.S.A. optimization digital 
computer program ISAFT. 



ii6 


/JOB A257 f'CPKI^^.LI^^S=b^i 

C 

C program tSAFF 

C DETcftfMAES CFTIFAL ^,C^STA^T GAIK fUTPLT ftECfAC< CtJ < I !OLL EP S 

C FCR FIMTi FIFE STATE RtGOLATCR RRCeiEFS 

C OLTPLTS ARc ASSEmEO TC K£ FIRST NF STATES 

C 

1 COMMON F(3,l),M(3,3),RC(3).tAh4I(},T),T,Kn,3),G(3,3l,^i\n.»,T), 

1 FT(3,31,FFTI(3,3) ,OF r(3,3,3,3),FCFTI(i,5,3,3l, 

2 C2FU3,3,3»3,3,3),DFi;FTI|3,3,3,3,2,3),FC2FI H 3 , 3 , i , S , 3 , i ) , 

3 VMS ) , VGRAC(9t9) ,GRACI (9,91 , S ( 3 , 3 ) f F 0 , 3 ) . R I 3 . 3 I ,K'S , \C f KF 

2 CIMENSICN A(3,3I,C(3,3) 

3 CCMPLtX*l6 F ,M ,RC 

A CnUBLE PRECISICK AF A T , T , K , G ■ Ch A I , F T , FF T I , CF T , F CF f 1 1 1. .'F I . CF OF II . 

1 FC2FTI ,VV>,VGKAL'tGRAOI 

5 COLBLc PRECISION GF ,CF2,TES 1 1 TES 1 1 .GKSTCF 

6 CCueUE PRECISION CFA*l,Dflt!S 

7 lu FORMAT lAFH.n 

8 11 FORMAT (ADL’.n 

9 IS FCRFAT (711.) 

10 2t FORMAT (lAIS) 

11 25 format ( • 1 ' , T9, ' ST a FES ' »AX, 'CCNTRCLS' ,AX > ' FEEOHACKS ’ , ax. ' I 3A IKS' ) 

12 30 FORMAT I//T2, • inverse GRADIENT MATRIX*) 

13 35 FORMAT ( // T 3 , * SYS T EF MATRIX A') 

lA 40 FORMAT (//IS.’NEl. CAINS') 

15 45 FORMAT t / / T3 , 'CCiN TRCL FAFRIX B') 

16 5C FORMAT {//I3.'GAt.\ MATRIX*) 

17 55 format I//T3,*TERFINAL TIFE T =*,F1H.7) 

18 60 FORMAT (//T7.'GA!N TCLERANCE ACHIEVED =',F13.7) 

19 62 FORMAT 1//T3, ’AVERAGE COST =',CZ .7) 

20 65 FORMAT I / / T 3 , * T =R F I M A L CCST FATRIX F'l 

21 7C format I//T3, 'RcCUREI) STCFPING TOLERANCE =*,F1''.7) 

22 75 FCRMAT I//I3, ’STATE WEIGHTING FATRIX (. ' ) 

23 8C FORMAT I / / T 3 , • CCN TRC L WEIGHTING MATRIX R*) 

24 85 FORMAT I * 1 • , ’ I T ER AT I CN MjFBER’,11-) 

25 90 FORMAT I / /T3 , ' SCLUT I CN IS CCFPLETE. A6CVE GAINS ARC CPTImAL') 

.26 92 FORMAT (//T3,’STEP SIZE TOO UAROE. NEW AVERAGE. COST WC.ULC HAVE PEL 

IN*,D16.7) 

27 93 FORMAT IT3,'CAIN ACJLSTFENT IS H.ALVrC) 

28 95 FCRFAT I / / I 3 , ’CON V£ RGENC £ TEC SLEW. preGRAM T HR F I N at ‘•-C ‘ I 

29 ICC READ (l,2';i N S , NC , NF , I G A I NS 

30 READ a, IE) I (A(I,J) ,J = 1,NS).I = 1.NS) 

31 ftpAo (1,13) ( ten , j) .J=i,NC) , i=i,ns ) 

32 READ (1,11) T 

33 READ (1,10 ( (F (1 , J) , J=1,NS) , [ = 1,nS ) 

34 READ (1,10 ( (G ( I , J) , J = l,NS ) , I = 1 ,Sa I 

35 READ (1,1. ) ( (R ( ! , J) ,J = 1 ,\C) , I = 1,NC) 

36 READ (1,11) GNSTCP 

37 WRITE (3,25) 

38 WRITE (3,15) NS , NC ,NF , I G A I N S 

39 WRITE (3,3>) 

4.) WRITE (3,1 ) ((A(I,J),J=1,NS),I = 1,NS) 

41' WRITE (3,45) 

42 WRITE (3,10 ( (B(i ,J) T J=1,NC) , I=X,NS) 

43 WRITE (3,55) T 
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h91T£ (3,65) 

‘r> 


k><IT5 (3,1,) ( (F ( I , J) , J^l.NS) , t-L,NSI 

A6 


)»R1TE (3,75) 

AT 


)>RITE 13, i;) ((Q(I,JI,J=1,NS),1 = 1,\S) 

AS 


V.R1TE (3,8:) 

A<i 


WRITE (3,10 t(R(I,J),J=il,^CI,l = l,NC) 

5y 


NFP1=NF*1 

51 


^FC = ^F*^•C 

52 


cr i:c3 1 = 1, ^s 

53 


CC ICCJ J=1,KC 

54 


Gti,j)=':.ccc 

55 

ICCC 

K( I,J)=0.CCC 

56 


IFdGAINS) 1203,12Cu,llCC 

57 

lli-C 

read (1,111 ( (< ( I , J) ,J = 1,AC) , I = l,f«F ) 

58 

12CC 

CCNTINUE 

59 


CG 1223 1 = 1, AS 

60 


CC 1223 J=5SF,AS 

61 


QHATd, J)=Cn,J) 

62 

122C 

AHATd,J)=Ad ,J) 

63 


1T = 1 • 

64 


WRITE (3,851 IT 

65 


WRITE (3,50 

66 


WRITE (3,i;) t(Kd,J»,J=l,.SC),I = l,NF) 

67 

1233 

CCNTINU6 

68 


cr 1250 1 = 1, AS 

69 


CO 125J J=1,AF 

7j 


A)-AT( I,J)=Ad ,J) 

71 


CHAT! I,J)=C( I ,J) 

72 


CO 1250 M = 1,aC 

73 


AHATd,J)=i6FATd,J»-Pd,M)*K(J,Ml 

74 


CC 1250 N2=1,AC 

75 

1250 

CHAT( I, J)=tHiT ( I,JI + Kd,M)«R(M,N21«K(J,A2) 

76 


CALL STRAP 

77 

1260 

call FEF.', 

78 


IT=IT*1 

79 


CALL GAINFMGF) 

80 


WRITE (3,62) GF 

81 


CALL MESCON 

82 


CALL GRAOAT 

83 


call Invert (vg8ad,gradi, AFC) 

84 


WRITE 13,3'.) 

85 


WRITE (3,10 dGRACId,J),J=l,NFC),I = l,AFC) 

86 


CALL N5)»RIT 

87 


WRITE (3,85) IT 

88 


WRITE (3,43) 

89 


WRITE (3,10 ((Gd,JJ,J = l,NC),I = l,NFI 

90 


IT1*C 

91 

1280 

CrNTINUE 

92 


DO 1300 [=1,AS 

93 


DC 13C0 J=1,AF 

94 


AHAT( I , J) = A d, J I 

95 


CHATd,J)=Cd,J) 

96 


DO 1300 M = 1,AC 

97 


AHAT(I,J)=AFATd,J)-8d,M)*G(J,til) 

96 


CO 13C0 N2=1,AC 



99 

13^9 

Qf^AH 1, J> =(.HA H r , Jl+cn ,M )»P 1 »0IJ,N2 1 

lOj 


C4LL STRAN 

101 


CALL GAIfv2t6»^2) 

102 


IF(Gf-~CF2,9T...'jC'’ ) GC TC 13- j 

103 


WRITE (3,92) GF2 

10<i 


WRITE (3,93) 

105 


CC 135) 1=1, ^F 

1C6 


CO 135’ J=1,AC 

107 

135v, 

G(I,J)=C( I,J)/2. + K(I ,J)/2. 

100 


WRITE (3, a:) 

109 


WRITE (3,10) l(G(I,J),J=il,NC),I = l,AF) 

lU 


ITl-ITl+l 

111 


IF(ITl.Lfc.5) GC TC 128^ 

112 


GO TC 3C0''. 

113 

138'o 

TESTl = ''.0n . 

ll<f 


CC 2CC) 1 = 1, AF 

115 


CC 2CC.3 J=l,Al- 

116 


IF(G( UJl.tC.^.CCC) GC TC lAt-j 

117 


TEST»CAeS(<G( I,J)-K(I,J))/G( I,J)) 

118 


r£STl*OFAXl( TuST.TESTl) 

119 


GC rc 2:c: 

12? 

19C- 

rF(G(I,J)-K(I,J).EC.(,.;CG) GC TC 2',:? 

121 


TESTl*lCCG.«GASTCP 

122 

200 C 

K(1,J)=G( I,J) 

123 


WRITE (3,6.;) IcSTl 

129 


WRITE (3,7,) GNSTCF 

125 


IF(TEST1-CNSTLF.GT.C.^-C;i GC TO 1?6C 

126 


WRITE ( 3,9 : ) 

127 


WRITE (3,62) GF2 

128 


GC TC 0CCC 

129 

3000 

WRITE (3,95) 

13J 

80 C w 

CCMTINUb' 

131 

^CCO 

CCNTihUE \ 

132 

133 


i-' 



119 


13-V 


135 


136 

■ 137 

138 

139 

I'iU 

1^1 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 
169 
17 J 

171 

172 
17.3 

174 

175 

176 

177 


SUBRCLTINE STR4f< 

C 

C COMPLIES IKE STATE TR4NSIT1CN H4TRIX 

C . _ 

COKMCN Ml 3, 3 I ,M(3,3J,RC13),AHAn3,3),T,K(3,3),G(3f 3)iUMAT(3,3), 

1 FTl3,3),FFTIl3,3UDFT(3,3,3,3),FDFTIl3,3,3»3»f 

2 D2FT(3,3,3»3,3,3l,BFDFTIl3,3,3,3t3t3),FCiFTI(3,3,3,Jj3,3)i , 

3 \(U19) ,'VGRAD(9,9) ,GRADI (9,9 I ,B( 3,3) , M 3 , 3 ) , R 1 3, 3 ) ,NS ,NC , \F 
CIMENSION «AAA19),RR(3) , R I (3 > , 4SCR (3 , 3 ) , A SO 2 I 3 , 3 ) , XR ( 3 I , X I ( 3 ) , VR I 3 

I I.VI OI .IANA(3).IRCW13.2I.VRN(3). VlN(3) ,Vl(3.41 

C0HPLEX*16 K,MI,RC 
C0HPLEX»16 OCMPLX.CCCNJG 

DOUBLE PRECISION AHAT , T , K , G , CHAT , F T , F FT 1 1 CFT , FC FT 1 , 0?FT ,„DF0FJ I, 

1 FC2FTI ,VW,VGRAD,GRAOI 

DOUBLE PRECISION 4 44 A , RR , R I , A SOR , ASC2 , XR , X I , VR , V I , VRN , V I N , R , V ECKGR 

1 ,VECMGI fVECMGSiSWi _ 

10 FORMAT 12D18.7I 

30 FORMAT ( / / T3 , • S VST EM EIGENVALUES*) 

CALL VECT(AHAT,AAAA,NS1 _ _ 

CALL HSBG(NS,AAAA,NS) 

CALL ATEIG(NS,AAAA,RR,RI,1AN4,NS) 

hRITE (3,31) _ . . 

WRITE (3,10 (RR{ I ) ,RI ( I ) , I = lfNS') 
call MSC(AAAA,NS,ASQR) 

DO ICO 1 = 1, N’S 

RCm=DCMPLX(RR(I) ,Rim I 
DO ICO J=1,NS 

ICO ASU2 ( 1 , J) nAStrt ( I , J) _ 

CALL EIGVECO ,AHAT,ASCR,W, IREW,XR,}('i,VR,VI,r'r( 1) ,RI(1),NS,NS,0, 

1 SW1,ITER,DIF,2) 

VECMGR = C.OCi; 

VEC»'GI=' .'.'C. 

CO lie 1 = 1, NS 

VECMGR=VcCMGRtVR ( I )»XR ( I J-VI ( I) «XI ( I 1 ^ 

lie VECMGI=V£CMG1+VR(I 1»XI U )+VI ( n«XRl I) 
VECHGS=VECMGR«VECKGR+VECMGI*VECMGI 

DC 120 1=1, NS . _ . 

VRN( I ) = (VR( I )«V£CPGR+Vl ( 1 1 *VECMGI) /VECPGS 
VIN(I) = (VI(I MV£CMGR-V,T( I } •VECHGI ) /VECHGS 

M ( I , 1 )=CCMPLX (XR ( I) , XI ( n } 

12U Mid, I I =DCMPLX( VRMI ) ,VIM I) ) 

CO ICC) R0LN’T = 2,NS 
K0UNT1=K0LNT-1 

IF(RR(KCLM)-RR (KCLNTl) I 2-..! , 14u, 2 j ) 

I4j CFNTINUE 

CO 15. [=1,\S 

M I I,KCUNT ) =CCLNJG (M I .KCUNTl ) ) 

150 MI (KCUNT, I l=DCCNJG(M.I (KCU-N'Tl ,11) 

GO TC ices 
2:0 CO 21. 1 = 1, NS 
CD 21J J = 1 ,NS 
212 ASO'U I , J ) =4SC2 ( ( , J ) 

CALL cIGVECl 3 ,AHA T ,ASOR , h , 1 RCW , XR, X I, VR, V I , RRIKOUNT ) ,RI (KfJUNT ) ,NS, 
1 NS,0,SH,ITER,CIF,2) 

vECMGR= ;c ; 



1,0 


17a V6CHGI=C.CC , 

179 DO 22C 1 = 1, ,NS 

180 VECMOR = VECr'GR+VR( n»XR t Il-Vl ( I) ‘XK I) 

181 220 VECMGl = veCr'GI*VR(n*XI ( I)+VI(I)«XR( n 

182 ^ECHGS = VECRGR»VECf'GR*VECHGI»veCt‘GI 

183 DC 230 r=l,\S 

184 VRN( I) = [VR{ T )»VECHGR*VI ( I I *VECHGI ) /VECHGS 

ies viN(i) = (viin •vechgr-vru (•vecMcn/vEcrGs 

186 H(i,KOUNT)=cCf'PLx[xRm ,xnu) 

_187_ _230 HI IKOUNT, n=DCMPLX(VRN( n,VIN( I M 

188 lOCb CONTINUE 

189 RETURN 

190 END 



191 


19?. 


193 

19^ 

195 

156 

197 

150 

159 

200 

201 

202 

203 

204 

205 

206 

207 

208 
209 
2l.') 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 


SUPRCLTlNc FEF^ 

C 

C COMPUreS RuCOIRED FLJ^CTIC^S oh the ST/!T-: TR415>II ICN ('Ar'UX 

c 

CGHf-CN P(3,3l,f'Il3,3),RCl3),AhAT(3,3),T,K(3,1),G(3,3),l.HATlJ,3l, 

1 FT( 3 , 31 *FFTI( 3 , 3 ),OFrni 3 , 3 , 3 liFCFTI( 3 . 3 i 3 t 3 ), 

2 n;2FT(3,3,3,3,3,3) ,CFCFTI (3,3,3,3,2,31 ,FC2FTI[3,3,3, 3,3,3), 

3 VV, ( 3 1 ,VGRAD(9,9I ,GkflD 1 15,5) , B( 3,3 I ,F ( 3, 3 ) ,R( 3, 3) ,\b,NC,t;r 
CIPENSICN CllJ, 3, 3, 3, 3, 31, 02(3,3,3, 3,3,3 ) 

C0HPLcX*16 tXl(31,ex2(3,3),Rl(3,31,R2(3,2), CDi'l , OOP 2, [.UP },UtP4, 

1 CLP5,CUP6,(3UP7 

.CnKPL£X«16 P,M ,RC 
CCMPLEX«16 CCcXP 

DOUBLE PRECISICN A[-AT,T,K,C,CHAI,Fr,FFTI,[:Fr,Fi:FTl,0?f F.uH.FT I, 

1 FC2FTI ,VV.,VCRAL,G.<ACI 

CCU3LE PRECISICN D1,C2 
NS1=NS-1 
NF1=NF-1 
NC1=AC-1 
DO ig; 1=1, as 
Exi ( I)=CCEXP (kc ( I ) «r ) 

CC IGl J=1,NS 

EX2( I, J)=CCEXF( (RC (I I+RCIJ) ) »T 1 
Rl(I,J)=RC{I)iRC(J) 

ICJ R2(I,J)=RCH)-RC(J) 

CC 1CC3 1=1, \S 
CG ICC 3 J = 1,AS 
FTn,j)=G.:c: 

FFTKI, J) = :.3C3 
00 SCO K1=1,AF 
CO SCO K2=1,AC 
CFT( I, J,Kl,K?) = C.GC'J 
FDFTI ( I , J ,K1 ,K2 )=0. :C'J 
DC 3C0 K3=1,AF 
DO SCO K4=1,NC 
D2FT ( I T J,K1 ,K2,K3 ,K4)=.5.C'D3 
CFDFTI ( [, J ,K1 ,K2,K3,K4)=;.:c: 

300 FD2FTI ( I,J,K1,K2,K3,K4)=G. D. 

CC ICGO N1=1,AS 

FT(I,J)=FT(I,J)+f'(I,M)»M(M,Jl«EXl(M) 

DO ICLO N2=1,AS 
CC ICCO N3=1,AS 

CUH1 = P(I,M)«PI(M,N2)«MJ,N3) 

FFTI(I,J)=FFTI( I , J l+CUPl«PnA3,(v2) "(EXPIM.ASI-I.I/kK Jl, )3) 

C'' 6C: K1 = 1,AF 
CC 6C0 K2=1,AC 

CLH2=P( I,M)*M (M,M2)«8(N2,K2)«P(K1,N3) 

IF1N1-N3) 320,310,320 

3 1C CFT( 1 , J,K1 ,K2 )=DFT ( I , J,Kl,K2 )-DUM2»P! (N3,J)*T*£X1(N1) 

GO TC 33 C 

320 QFT(r,J,Kl,K2) = DFT(I , J , K1 ,K2 1 -DUP 2 *P I ( .S3 , J ) • ( UX 1 ( N 3 1 - EX 1( 'i 1 1 ) /R 2 ( M 
13, Nil 

33C CO 6C0 N4=1,SS 
CO 6C0 N5=1,S' 

DUM3=0UPl»PHS3,S4)*B(N‘4,K2)«P(Kl,Ai51»PI(S5,N2) 




237 

238 

239 
2^0 

241 

242 

243 

244 

245 

246 

247 

248 

249 
2 5C 

251 

252 

253 

254 

255 

256 

257 

258 


259 

263 

261 

262 

263 

264 

265 

266 


267 

268 


269 

27*J 

271 


272 

273 


274 

275 

276 


IF(N3-N5) 35;,34:,35J 

340 FCFTI([,J,K1,K2)=FCFTI([ , J , K 1 , K 2 ( -lMJt'3 « ( E X2 ( M , i\3 1 » ( k U a,N3)»T-l, 

1) +l.)/tHl(M,N3)»Rl(M,63n 
Gf TC 36U 

35C FDFT t ( I ,J,K1,K2 )=FCFT1 ( ! , J tK 1 , K2 l-t!UR3* ( « EX2(Nl,N5)-l.)/r<liNl,\5)- 
■1(EX21N1,N3)-I.I/R1(N1,N3J »/K2{K5t N3) 

360 CO 6CC K3*1,KF 
CO 6C0 K4=1,NC 

DtH4 = 2. •nUF2«M (^3 ,N4)»H(N4,K4)»F(K3,N5)*M(N5,J) 

IF(N3-N5> 4-;:,370»4..f 
37C IF(N1-N3) 39.;,380.39C 

380 C2FTt I , J,Kl,K2,K3tK4)=C2rT(I , J,K1,K.2,K3,K4I+CUF4«T»T«EX1(N1)/2. 

GC rc 455 

39C C2FT(I,J,KltK2»K3.K4)-C2FT(IfJ,Kl,K2,K3,K41+CIJ>'4»t(R2{.x3,-'ill«T-l.) 

1«FX1(N3 1 + EX1(M))/(R2(N3,M)»k2IN3,M)) 

GC TC 455 
4C0 CCNTiNUE 

IF(M-N3) 42:,41C»42C 

410 C2FTI I , J,K1,K2,K3,K41 = C2FT( [ , J , K 1 , K2 , K3 , K4 ) -CUF4»T«EXH M )/R2 ( ^5 
111 

CO TC 43C 

420 C2FT(I,J,Kl,Kc,K3.K4(=r2FTU,J,Kl»t<2TK3,X4)-CUF4«(EXl(N3J-6Xl(Mn 
l/(R2tN3,Nll•Rx(^5.^3)l 
43C IF(N1-N51 44:,45C.44r 

443 C2FT ( I , J ,K1 f K2,K3»K4 l=C2FT< I , J,K1, K?,K3 tK4 I +[0^4 » ( EX H N 5 )-cX 1 ( N 1 1 ) 
l/(R^(^15,M)•t^^(^5l^3U 
Gr TC 455 

45C C2FT{I,J,Kl,K2,K3fK4)=C2FT( I , J , M, K2, K3 , K4 1 *CUR4*T • EX 1 ( .n 5 I /R2 ( \S,fv 
131 

455 CO 630 N6»1,N5 
CO 6C3 N7=1,NS 

DUH5 = RI tN4,^5 1•e(^5,K4)*^'(K3,^6)*^'I(^6,^7) 

CU’46 = CUF2*F I (N3.N71*F'(J,K4)»CUF5 

CUH7 = 2.«Ff I,^l)•M(N1^^7)»^(J,^21*^'I(^2,^^l•el^3,K21«^'(Kl,N4) •CLR5 
IF(Ni-N3) 49C,460»49C 
460 1F(N4-N6I 48j,47C»48C 

470 DFDFTK I , J t K1 ,K2 , K3 , K4 1 =CFCFT 1 1 I , J ,K1 , K2 ,K3,K4 > ♦CUF6* t ( T»T»RHN1,N 
14)*Rl(Vl,\4)-2.»T»Rl(M,N4) + Z,l«i;X2(M,K4)-2.)/(Rl(M,'J4 1»Tl(Nl,N4 

2) *R1(N1.N4) 1 
GC TC 52C 

480 CFDFTI t I , J , K 1 ,K2 , K 3 , K4 ) =CFCFT I t I . J t K 1 , K2 . K 3 ,K4 )+CUf'6«( t t Bl(Nl,N6 I * 
lT-1, )*EX2(N1,N6) + 1. I/(R11N1,N6) •Rl(^!l,^6} )-l IR1(M,N4)*T-U)*EX2(N 
21,N4) + 1, ) / (RHM»N4) «RHM,N4 I ))/R2tN6,K4 ) 

GO TC 52o 

496 IF(N4-N6) 51':,50C»51C 

5CC CFDFT I ( 1 , J ,K1,K2,K3 ,K4)=CFCFTI I I , J , K 1 , K2 , K3 ,K4 ) +0UW6» U (^<l^^3,N4 I • 
lT-1. )*EX2( N3,N4 )+l. )/(Rl (N3,N4) «R1 (N3,N4 I )-( (Rl(Nl,hi41*T-l. )*eX2(N 
2l,N4) + l.)/tRi(M,N4)«RltM,N4n )/R2(N3,M> 

C-0 TC 520 

510 CFDFTni,J,Kl,K2,K2,K41=CFCFTHI,JtKl,K2,K3,K4>*DUf'6«((ex2(.N3,.N6l- 
ll.)/R1^^3,^6l-[EX2(N3,^41-l.)/Rl{N^,^4)“(EX2tNl,N6|-l.)/Rl(Nl,N6l+ 
2{ EX2(N1,N4)-1 . l/RUMaAH / l!62(N3,»ll) •R2(\-6,N4 I ) 

520 IF(N4-N6) 56Ci530f56ij 

53C IF(N2-N4) 55'.,54C»55G 

540 F02FT I ( t . J .Kl ,K2 , K3 , K4 I =;FD2FT I ( I » J i K1 f K2 , K3 ,K4 I + , 5»CCH 7 » I EX 2 ( M ,N 2 



J>3 


l)«(T«T»RllM,N2l«Rl(M,\2)-2.«r»«l(Nl,N2) + 2.)-2.l/(‘^UMtN2)»Kl(M 
2,N2) »R1(M ,N2) ) 

277 GO TC 6C: 

278 5S0 FD2FT I tl , J,KX jK 2 , K3 , ) =FC2FT I t 1 , J , K I , K2 , K3 , > ♦CUI'V* ( 1 F?(2(M tN'. I • 

l(Rl(M,N4)«T-l.) + l.)/tRHM,SAI»RllM,N^n-[EX2{Nl,N^)-l.}/{PHNl, 
'2N4)*R2(N4,f<2> ) + (EXZ(M,N2)-l. )/(Bl(M,N2)»R2(NA,N2) ^/ft2<NA,^2) 

279 GO TC 6CO 

28U 56C IF(N2-NA) S8'.,S70,58C 

281 57o FD2FTK I , J , K I ,K2 , K3 , KA ) =FC2r T i ( I , J , K 1 . K2 ,K3 , KA ) -DUK 7 • ( ~ X 2 ( M ,N2 ) • ( 

lRHNl,N2)»T-l.) + l.)/(RltM,K2)»«llMtN2)*R2IN6,N2)) 

282 GO TC 590 

283 580 F02FTI ( i,J,Kl,K2,K3,KA)=;FC2FTItI,J,Kl,K2,K3,KA 1-CUR7«< I tX2(M,NA) 

l-l.)/RUNl,NA)-t£X2(M,N'2)-l. )/RllNl,N2) )/{R2(NA,K2 1*R21Nfc,MA 1) 

28A 59C IF1N2-N6} 59A,592.S9A 

285 59 2 FOZFT I( I, J,K1 ,K2 , K3 , KA ) =F D2FT tl I , J , K I , K 2 ,K3 ,K A I +DUR7* 1 1; X 2 ( ,M I , Nit ) • 

l(RHM,N6)»T-i.) + l. ) / (Ri(Nl,K61 •R2^^6,^A)•RHNl,^!0| ) 

286 GO TO 6C0 

287 59A FD2FTK I , J ,K1 ,K2 tK3 , KA 1 =FC2FT I t I , J , K 1 , K2 .K 3 , KA I +0UK7»( ( £X2tM,Nfc)- 

U, ) /Rl(MiN61-(EX2(Nl,N2l-l, I/R1(M,K2) ) /1R2(M6,N2 l*R2(iN£,i<A ) ) 

288 6Cti CONTINUE 

289 CO 7C0 Ki=l,NF 

290 CC 7C0 K2=1,NC 

291 DC ICC K3^1,KF 

292 on 7C0 KAaI.NC 

293 Dl(I,J,Kl,K2,K3,KA ) = C2FT( I , J , Kl , K2 , K 3 , K A ) 

29A 700 D2( I,J,K1,K2,K3,KA)=FD2FTI tt,J,Kl,K2,K3,KA) 

295 on 8C3 Kl=l,NF. 

296 DO SCO K2=1,NC 

297 DO 800 K3=l,NF 

298 on 800 KA=1»NC 

299 D2FTt I, J,Kl,K2fK3,KA)=(01( I.JtKl,K2tK3,KA)+CUI,J,K3,KA,Kl,K2)l/2. 

300 800 F02FTI (I,J,Kl,K2,K3,KA)=(D2{t,J,Kl,K2,K3,KA)+C2(I*J,K3,KA,Kl,K21)/ 

12. 

301 ICCC CONTINUE 

302 RETURN 

363 END 



124 


304 


305 


3C6 

307 

308 

309 
313 

311 

312 

313 

314 

315 

316 

317 


SfQHCLTINE ■J*1^F^(CF) 

C 

C COF'PUTE.'i iv£«4G£ CCST 

C 

COPFCN FO,’) fPI (3,31,RC(3) ,«Hfln3,3l,T,K(3,3 J,G(3,3),(.>iiTt3, 3) , 

1 FT(3,3 },FFTI(3,3) ,DFT(3,3,3,3} ,FCFTK3, 3,3t 5). 

2 02FT(>,3,3,3,3,3) ,CFCFrl.(3,3,3,3,3,3) ,FC2FTn3,3,3, J.'s, i) , 

3 VK(9) , vGR/SC(9,9 1 ,Gft AO I < 9 , 5 1 , P ( 3 , 3 ) , F ( 3 , 3 ) , rs ( 3 , 3 ) , S , '.C , NF 
CnHPLex»l6 F,M,RC 

D0U3LE PRECISICN AF AT , r i K , G , CHAT , F f , F F T I , CFT , F{, FT I , GIF T f t FOT fl, 

1 FC2FTI ,Vh,VGRAt:,GRAC! 

DOUBLE PRECISICN GF 
G F “C • CD Z 
CO ICCG M = 1,AS 
DO ICCG N2=ltAS 
GF=GF+0FAT(M,S'2)*FFTI (N2»N1) 

DO ICOT ^3=1,^S 
icoa GF=GF+F<N1,A2)*FT(A2,N3) 

GF=GF/NS 
RETURN 
END 




318 


SUBRCUTtNE ^esCO^ 


319 


320 

3.21. 

322 

323 
32A 

325 

326 

327 

328 

329 
3 3v> 

331 

332 

333 
33A 
135 

336 

337 

338 

339 


C 

C CCHPUT6S NtCESSiRY CCNDITICNS 

c. 

CnMHCN H(3,3) ,NI(3,3),8Ct3),4hATt3,3),T,K(3,3J,GI3,3),iJhATI3,3l, 

1 FT(3,3),FFTI(3,3),DFl(3,3,3,3)fFCFTn3,3,3,3), 

2 D2F713,3,3»3,3,3),DFDFTI(3,3,3,3,3,3),FC2FTn3,3,3, J,3,3>, 

3 VW(9 ) ,YGRAD (9,9I.,6RACI (9,9|,B(3,3),F(3,3>,R<3,3I,NS*NC,KF 
COMPLEX-16 K,MI,RC 

OQUBLE PRECISXCK AHAT , T , K , GjCHAT , FT , FFT I , CFT, FCFT [ , D2 FI , CFDF T I t 

1 FD2FTI,VH,VGRAD,GRA0I 

10 FORMAT (//T1, 'NECESSARY CCNOITICNS VECTOR' I 
20 . FORMAT (ADia.7) 

NFC=NF*NC 
DO 1103 1=1, NF 
DO IlOO J=il,NC 
IN=NF«t J-l»tI 
VM(IM=C.003 
00 1100 Nl=l,hS 
CC icco N2=1,NS 

VW( IN)=VhnM + 0HAT<Nl,N2)«FCFT! (N2,M, I . J) 

DO ICCO N3=1,NS 

■' 1000 VViflN) = VHt iM*FTAl,N2)»FT(N2,N3)«CFT(Nl,A3,I,J) 

CO IlOO N4=1,NC 

_ 1100 VhnNJ = VWnN)4R(J,N4)*K(M,N<t)*FFTI(Mjll 
ViRifE' (3,10) 

hRITS (3,2:)' (VKU) ,I = l,NFCJ 

RETURN _ . 

END ■ ' 



340 


SLBRCtTINE GRADM 


I. 6 


341 


342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 


361 

362 

363 

364 


C 

C COMPUTES GK4DIEM MATRIX 

C 

COKMCN I3,3» ,MI (3,31 ,BC(3) ,AMftT{3t 3) ,T.K(3,31 ,G( 3t3),t.HAT(3,3 > , 

1 FT(3,3>,FFTI(3,3),DFT(3,3f3i3) ,FCPTn3,3i3,3J, 

2 D2FT(3,3,3»3,3,3),DFDFTI(3,3,3,!i,3,3),EC2FT!(3,3,3,3,3,3), 

3 VhI9) t vGRA0(9,9l ,GRA0M9,91 , Rt 3 , 3 ) i F ( 3, 3 ) , « ( 3, ? I ,^S .«C ,.\t- 
CCMPLEX»16 F,MI,BC 

DOUBLE PRECISION AFAT , T , K , G , CM AT , FT , FF T I , CFT ,F OF F 1 , 02FT , CF OF J I , 

1 FP2FT I tVWTVGRADtGRACI 

10 FORMAT (//T3, 'GRADIENT matrix*) 

20 FORMAT (4D13.n 
NFC=NF*AC 
CO 11C3 I=!l,NF 
DC IICJ J=il,NC 
IN»NF*I J-1 )♦! 

DO 1100 K1=1,NF 
DC 11C3 K2*lfNC 
I0=NF*tX2-l)*Kl 

VGRACI IN, ID)=R(J,K2) *FFTI (K1,I» 

CO IICO NI*1,nS 
CC lOOO N2=<l,NS 

VGRADnh,IC)=FVGRAD(IN,ID)+CtIATtNl,N2<*<CFCFTt(N2,N|l,Kl,K2, I,J)+FD2 
IFTHN2,M,K1,K2, I , J) » 

DO i:Cj N3=1.NS 

1000 VGRADI IN, lD)f VGRAC( IN,ID)+F(M,N2)»(DFTIN2,M,Kl,K2)fOFr<Vl,N3,I,J 
2) ♦FT{N2,N3 )402FT{M,N3,!<1,K2,I , J») 

00 llOO N4=1,NC 

IICC VGRACdN, IC) = \(GRAC[IN,ID)+R<K2,N4) (N 1 , N4 ) • I FE FT 1 1 N1,K 1, I , J ) ♦ 

1 FCFTI (K1,M,I,J) ) ♦ RU,N4)«K(N1,N4)? 

2 (FDFTIIN1,I,K1,K2)+FDFTH I,M,K1,K2)) 

WRITS (3,10) 

WRITE ( 3 , 20 ) ( (VGRAD(I,JI,J=1,NPC) , I=1,NFC) 

RETURN 

END 



36 & 


366 

367 

368 

369 
37j 

371 

372 

373 
37A 

375 

376 

377 

378 

379 
38 3 

381 

382 

383 
38^ 

385 

386 
3E7 

388 

389 

390 

391 

392 

393 

396 

395 

396 

397 

398 

399 
6CG 
6C1 
6C2 

6C3 

606 

605 

606 
607 
6C8 


SOBRCUTlNe lNveR7{ft,P,M 
C 

C INVERTS A 1C GIVt 8 

C 

DIMENSION AI9,9I ,8(9,9) 

CCUBLE PRECISION A,e,C,0,X 
COUBLE PRECISION DABS 
IF(N-l) 1CJ,1:-J,1C1 
ICO 8(1, 11=1. coi;/Ati,i) 

RETURN. 

ICl CO 1C2 1=1, N 
DO 1C2 J=1,N 

102 8(1,31=0.000 
DO 1C3 1 = 1, N- 

10 3 8(1, 11=1, :cc 

C PICK UP PIV.CT ELEMENT 

CO 116 K=i,N 
L=K 

IF(N-K) li:,U0,lC6 
106 I=Ktl 

CG l':6 JJ = t,N 

IF(DABS( A( JJjKl )-DA8S(A(L,Kl mC6,lC6., i05 
105 L=JJ 
1C6 CONTINUE 

IF(L-K) 107,110,107 
C PERFCRH RON INTERCHANGE 

1C7 CG 108 J=K,N 

C=A(K,J1 , . . . _ _ 

A(K,3)=A(L,jl 
lOa A(L,J)=C 

CC 1C9 J=1,N 
C=B(K, J 1 

e(K,Jl=8(L,Jl 

109 B(L,J1=C ■ ■ . _ 

C COLUMN ELIMINAflCN 

110 DC 116 1=1, N 
IF(K-n lll,li6,lll 

111 C=A( I,K1/A(K,K1 
CO 112 J=K,N 

112 A(I,J1=A{ I,J)-C*A[K,J) 

A(I,K)=C.CDC 

DO 113 3 = 1, N 

113 BII,J1=B(I,J)-C»B(K,J) 

116 CONTINUE 

C SOLVE FCR INVcRSE 

CC 115 3 = 1, N 
CC US 1=1, N 
X=Bn,J)/A(I,I) 

115 BtI,J)=X 

RETURN 

END 



4C9 


SUBRCUTINE NEnRH 


1. 8 


410 


411 

412 

413 

414 

415 

416 

417 

418 

419 
4ZJ^ 

421 

422 


C 

C PERFORMS Ncl)TC6 R4PHSCN ntRATICN 

C 

COMMON M(3,3» ,MI (3,3) tRCOl ,3HA7f 3f 31 , T, K (3, 3), G(3, 3 ), ■.'HA T( 3, 3), 

1 FT(3,3J,FFTI(3,3),CFT(3,3,3,31,FCf 71(3,3, 3,3), 

2 C2FT(3,3,3,3,3,3),0FOFTir3,3,3,3,3,3),FC2FT,I.(3, 3,3, 3,3,3 1, 

3 VR(9'),VGR(lC(9,9i,G«A0n9,9),e(3,3),F{3,3l,Rl 3 , 3 ) ,'N'S, NC , NF 
CCHPLeX«l6 H.MJtRC 

DOUBLE PRECISION AhAT, T,K ,G , CFiAT , FT , FF 7 [ , Cf T ,FLF f I, C2F T , CFDFT I, 

1 FC2FTI,VM,VGRAC,GHA0I 

CO ICOO I=il,NH 
DC icon J=1,NC 

G[1,J)=KU,J) 

IK=NF»(J-1 )♦! 

CO -1C33 K1 = 1,NF . 

DC ICCg K2«1,NC 
IC=NF*(K2-1 I*K1 

lOCC G(I,J)=GI I,J)-GRAOI(IN,ICJ»VK( rp) ' . 

RETURN 

END 



129 


423 


424 

425 

426 

427 

428 
42>9 

430 

431 


SU8ftCUTIN£ VECT(/th4T»4«AAfNSI 
C 

C CCNVERTS 4HAT TO SINGLE SU8SCRIPT fORN AAAA 

0 T M E N sT ON ' Tha'tTb Vi ) tA 4* iT('9i 

DOUBLE PRECISION AHATtAAAA 

00 ICO J=1jN5_ 

CO ICC 1=1, NS 
K = (J-H»NS+-I 

ICC AaAA(K)=AHAT(I ,J) 

RETURN 

END 




432 


SUBRCUTINE G4IN2(GF21 


' 433 


434 

435 

436 

437 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 

456 

457 

458 

459 

460 

461 

462 


C 

C COMPLIES AVERAGE COST FCR CCNVERGENCE CHECK 

C, . . 

C'OMHCK K(3,3> ,MI(i,3),RQ43» ,AHAT (3,3),T,K(3,3),GJ3,3I,0FAT(3,3 I, 

1 FT(3,3),FFTn3,3»,0FTJ3.3t3t3l , FOFTI ( 3 , 3 , 3t 3 J , 

2 n2FTt3t3.,3j3,3,3),'0FDFTU3t3,3,.3,3,3|,FC2FTI(3,3,3,3,2,3), 

3 Vh<9»,VGRAD(9,9) ,GRADI(9,9),ei3,3),F{3,3),R(5t3) .6S,NC,\F 
C0HPLEX*16 EXl(3),EX2n,3) ,R1 ( 3, 3 ) , CUMl 

C0J1PLEX*16 t'jPirAC 1 

CbMPLEX*l6 "COEXP 

DOUBLE PRECISION AHAT, T ,K,Gf CHAT ,FTt FFT 1 ,CFT , FCFT 1 1 C2FT , CFDFTI , 

1 FC2FTI,VH.VGRAD,GRACI 

DOUBLE PRECISION GF2 
DD 5CC 1 = 1, NS 
EXim=CDEXP(RC(I)«T) 

DC 5C0 J=1,NS 

£X2(I,J)=CDEXP( IRCm+RClJ))»TJ 
5CC Rll I,J)=RC(I)'»RC(J) 

CO ICOO 1=1, NS 
CO ICOD J^l,NS 
FT( r,JI=G.:0-3 
FFTI(I,J)=G,0C3 
DO ICOO N1=1,NS 

FT(I,J)=FT(I,J)+Hn,fil|«MI(,M,J)#EXllNl) 

CO ICO'J N2 = 1,NS 
DO ICCO N3=l,NS 

DUML = MI I (M,N2I*M J,N31 

ICOO FFTin,J)=FFTII I,J)+CUPl*'PnN3,N2)«lEX2(M,M3)-i,)/[U(Nl,N3) 
GF2=C.00G 
DD 2CC0 M = 1,NS 
DC 2CC3 N2=1,NS 

GF2 = GF2 + QHATIM-,N2I »FFTI (N2,M) 

CO 2CG9 N3=1,NS 

2000 GF2=GF2 + F(M,n21»FTIN2,N3)«RTIN1,N3I • 

GF2=GF2/NS 

RETURN 

END 



463 SUBROUTINE HSBG(N>A»IA) 

C 

C CONVERTS A TC UPPER H6SSENHEHG FCRf* 


C___. 


> ^ 

O 

0^ U1 > 

DIMENSION A(9) 

OOUBtE PRECISION A,PIV,T,S 
DOUBLE PRECISION DABS 



467 

L=N 



468 

NIA=L»IA 



469 

LIA=NIA-IA 




20 IF(L-3) 360.40f40 



471 

40 LIA»LIA-IA 



472 

L1=^L-1 



473 

L2=L1-1 



474 

ISUB=LIA+L 



475 

IPIV=[SL'B-IA 



476 

PlV>DA8S(A(IPlvn 



477 

IF(L-3) 90,90,50 



478 

50 H»IPIV-IA 



479 

DO 8C I=U,M,IA 



480 

T=0A8S(A( I 1 ) 



481 

IFIT-PIV) 80.80,60 



482 

60 IPIV*I 



483 

PIV»T 



484 

BO CONTINUE 



485 

90 IF(PIV) 100,320,100 



4 86 

100 IF(PIV-CABS(A(ISUB)n 180,180.12C 



487 

120 H=IPIV-L 



- 488 

DO 14(1 l3i,L 



489 

J«M+I 



490 

T=«A(JI 



491 

K=LIA+I 



492 

AU)=A(K) 



493 

140 A(K)=T 



494 

M«L2-H/IA 



495 

00 160 I>L1,N1A,IA 



496 

T.AIl) 



497 

J*I-M 



498 

A(I)=:AIJ) 



499 

160 A(J)=T 




IBO DO 200 I-L,LIA,1A 



SOI 

200 Am^ACn/AIISUB) 



SJ}2 

J«-IA 



503 

DO 240 1=^1, L2 



504 

J=J+IA 




505 LJ=L4J 






515 

516 

517 

518 

519 

520 

521 

522 

523 
52 A 

525 

526 

527 


.1 '2 


S=A{LK) 

LJ=L-IA 
CO 260 J=1,L2 

JK=K+J. 

LJ=LJ+IA 

280 S=S+A(LJ)*A(JK)*1.C0C 
300 A(LK)=S 

DO 31C I=L,LIA,IA 
310 A(I)=O.CDC 

320 L- Ll 

GO TC 20 " 

360 RETURN 
EN_D 



133 


SU8RQUTINE ATEIG tf'tA>RR».Rl*lANA»lAI.: 


C GOHPUTES ROOTS OF UPPER HESSENpERG. HATRI^ 


, DOUBLE PRECtSrON DABS, OSflRT, CHAX1\, - 

INTEGER p, PI, 0 

E7=i.QD~a --• "r . 


; 



■ MAXn=30 
N=H 

20 N1=N-1 



30 NP=N+J. 

IT=0 

DO 40 1=1,2 


549 P AN1=0.000 , ' “ v:-. 

550 R=c.000 - ' 

551 s=o.oDO • ' '■■ ' •■ 

555'' " N 1 N=IN+N 1 ; Y..V-' 

556 . 

557 60 ..- 











579 

580 

581 

582 

583 
5 84 

585 

586 

587 
583 

589 

590 

591 

592 

593 

594 

595 

596 

597 

598 

599 
6CJ 

601 

602 

6C3 

6C4 

605 

606 

607 

608 
609 
61 j 
611 
612 

613 

614 

615 

616 

617 

618 
619 
62. 
621 
622 

623 

624 

625 

626 
•627 

62« 

629 

63j 

631 


9R(N)=AtNlM) 

13C R1(N)=%CDC 
R I (M )=C.CCC 

RKNIUC.-. . 

GO TC 16C 
140 RR(N1)=L 

RR(N)=U . 

RI(N1»=V 
RI(N)=-V 

16C IF(M2U28C,12£'J,ieC , 

180 MN2=N1N1-!A 

RH0D = RR(N1) «RH(M)+RI tMURKNU 

EPS=E 1 J»DSCRT (R 8 CC) 

IF^DABS( 4 tN 1 ^ 2 ) )-EPS) 1280,1280,240 
240 IF(DABS( AUM) )-ElC*DARS( AtKM n 1320 , 13GP , 2 50 

25C IFt0ABS(PAM-AtMN2))-CAfiStA(NlN2) l?£6)._124C,_124Q,26C , 

260 IF(DABS(P4N-AtKM) )-CABS(A(NM) 1«E6» 1240,1^40,300 
3C0 IFtn-MAXITl 320,1240,1240 

320 J = 1 ' . ■ ... 

CO 36P t=l,2 
K=NP-I 

IF(OAB5(RRtK)-PPR 1 1 ) KRABSIRI (KI-PBU I) »-CELT Ax{ DA8S (RR (K 11 
1 +CABS tRI (K I 1 n 343,360,363 
340 J=J+I 
36C CCNTINU6 

GC TC {44C,'i6..,460,4SC» , J 
44C rt = 3.CC‘3 
S=O.CD9 
GC TC 50C 
46C J=N+2-J 

R=RR{ J) •RRU) 

3=RR(J)+RRt J) 

GC TC 5FJ 

48C R=RRtN)»RRtM)-RI tM*RI(M) 

S = fiR(N}-*RR(M> 

5CC PAN = A(NM) 

PAN 1 =A(NIN 2 I 
CO 523 1 = 1,2 
K=NP-I 

PRR( I )=RR (K ) 

52C PRI n I=R1 I ■<) 

P=N2 

IF(N- 3 ) 6 ':C , 62, ,525 
525 IP 1 =N 1 -I 2 

CG sec J=2,N2 
IPI=IPI-IA-1 
IF(DA8S(AnPI))-£PSl 6F J , 6C: , 53'" 

53C IPIP=IPI+IA 

IP IP2 = tP(P + TA 

C = A(IPlD(*tAUPIF)-S)FAnPIP2)tAnpiPfl-HR 
IF(0>54C,56' ,;4-r. 

54 U IF(DABSIA(IFU«A(IPIP-fl) ) * t CABS ( A II P I P ) + i 1 1 P I P 2 + 1 1 1 *C A? S( A < IP IP 2 

I + 2 )) 1 -C«BSIC)*HPS 1 62 ;, 6 Z.', 56 -■ 

56, P=fU-J 
58 J CCNTIN’JE 




632 

6CC 

e=p 

633 


GO TO 680 

63A 

62u 

P1=P-1 

635 


0=P1 

636 


IF(P1-1)68C,680.65C 

637 

6 50 

DO 660 I»2tPl 

638 


IPI=IP!-IA-1 

639 


IF(OABS(A( IPI) 1-EPSl 680»680,66C 

640 

660 

C=Q-1 

641 

630 

II=(P-1>*IA+P 

642 


CC 1220 I = P,M 

643 


Ml=n-IA 

644 


IIP^II+IA 

645 


IF(I-P)72Cf7Gu,720 

646 

7CC 

IPI=I1'H 

647 


IPIP=IIP+1 

648 


GlsAdl )»(A(II)-SUAniP)«A(IPI)-»R 

649 


G2=At IPH*( A(IPIP)+A(in-SJ 

650 


G3=A(IPn*AfIPIP + l| 

651 


A(IPI+l) =C.OCD 

652 


60 TC 783 

653 

72C 

Gl=Anil) 

654 


G2=Anil+l) 

655 


IF(I-N2)74C,740,76C 

656 

74 C 

G3=Atni + 2) 

657 


GO TO 780 

658 

760 

G3=O.ODC 

659 

78(t 

CAP=CSQRT(G1»G1+G2«G2*G3»G3) 

66U 


iF(CAp)aao,86C,80o 

661 

800 

IF(G1)82C,84Q,84Q 

662 

620 C4P=-C4P 

663 

E4C 

T-Gl+CAP 

664 


PSI1=G2/T 

665 


PSI2=G3/T 

666 


ALPHA = 2.0DO/'ll.OCO + PSU*PSIl*PSI2*PSt2) 

667 


GO TC 880 

668 

863 

ALPHA=2.0D3 

669 


PSI1 = 0.0DC! 

6 7-j 


PSI2=0,<JCC 

671 

88 :• 

IF(I-C)90O,96w,9C0 

672 

90 C 

IF(I-P)92C,540,920 

673 

92C 

A( ni)=-CAP 

674 


GO TC 960 

6 75 

941 

a ( II 1 1 =-A ( I (1) 

o76 

96: 

IJ = I1 

677 


CO l';43 J = I.N 

678 


T=PSI1*AI IJ+1) 

6 79 


IFn-M)983,lCQO,100C 

680 

980 

lP2J=IJ+2 

681 


T=T+PS12»A (IP2J) 

682 


ETA=ALPhA* ( T+al IJ ) » 

683 


A { IJ l=A( IJ )-cTA 

684 


A(1J4 1I=A{ IJ+U-PSIUETA 

685 


IFn-Nl)132';tl340,lu40 

686 

1C2C 

AUP2J) =A( IP2J1-FSI2«ETA 



687 1040 IJ=IJ+IA 

688 IFtI-Nl)10e!:,lG6C,lC60 

689 1C60 K»N 

690 . .. .GQ.TC. LlCjt. 

691 1080 K*I+2 

692 1100 IP=IIP-I 

_ 693 00 1180 J?.q,il( 

694 jiP=IP+J 

695 JI=JIP-I4 

696 _ T=PSI l»A(jrP> 

697 iF(f-Nl) 1120,1140,1140 

698 1120 JIP2=JIP+I* 

699 T=T+PSIJ*AUI_P2J ' 

70.) 1140 ETA=«ALPHA<*IT*a‘( jf) ) 

701 A(JI)=AtJn-EI4 

702 AtJlP)=A< J1P)-ETA»FSI1 

703 1F(I-N1U16C, 1180,1180 

704 1160 A( JIP2»»A(JIP2I-ETA*PSI2 

701 1180 CONTINUE 

706 IF{I-N2h20C, 1220, 1220 

707 1200 Jl=II+3 

708 „ JJP=JI+iA 

709 jlP2=jlP+IA 

710 ETA»ALPHA*PSI2«AU1P2) 

J.11. _i<J_n=-EXA 

712' A(JIP)=-ETA«PS1 1 

713 A( JIP2)=A(JIP2)-STA«PS12 

714 1220 11 = 1TP*1 I _ 

715 IT-'n + l 

716 GO TO 60 

7l_7 1240 lF(DABS(AlNNi n-DAeS(A(MN2M» 1300,1280,1230 

718 1280 IANAtN)*C 

719 IANA(NU=2 

720 N=N2 _ 

721 IF(N2J14C0,14Cd,20 

722 1300 RRtN)=AtNNJ 

^723 _RIIN)_50,JID,5 

724 IANA(Ni=l 

725 IF(N1)14C0, 1400,1320 

L26 132C N=N1 

727 GO fc 2C ” 

726 1400 RETURN 

729 END 



730 


SUBRCUTINE HSC< ^>)<>/>iNSyASQR) 
C 

C COMPLieS ASCR = 4iAA*AAAft 

c 


731 

732 ■ 

733 

DIKENSiON AAAAt9)yAS0M3.3} 
DOUBLE PRECISION AAAA.ASCR 
CO 100 lal.NS 

«r< 


734 

DO ICO J=l,NS 



735 

ASOR(ltJ>^C.CCO 



736 

DO ICO K=1,NS 



737 

Kl=NS*I J-ll+K 



738 

K2=NS«(K-ll*I 



739 

100 AS0R(I.J)qASQRiIyJ)«AAAAtKll«AAAA(K2l 


* 

740 

RETURN 



741 

END 





138 


742 


SUBRCLTIKE Eiovecnvc, i, 0, h, IRChr XR , XI. VR, VI, RCUfKE, 



1 

RCGTIE, NE, ^^'AX, T2 , SUl , CCOME, ERH.f'R^') 

i S' 


C 

SLBRCLTINS TC EINC THE EIGEMVECrnRS CF A NOK-SYRhETR IC MATRIX 

■: 


C 

BY A RGCIFIED VtlUKlNSCA'S INVERSE ITEFATICN RETI-LC. 

r:S‘ 


c 

CCNTROt IVC CCCE IS 

h S' 


C ■ 

1 FIND CNLY THE REGLLAR EIGENVECTERS (A X = LAKBOA X) 



c 

2 FINC CNLY the TPANSPCSEC EIGENVECTORS lAT V = LAMBDA 

V 1 tfS* 


c 

3 FIND eCTH TYRES CF E I GEN VEC TC P S . ' 

i 5‘ 

743 


CIMENSICN A(3,3),e(3,3),M3,2),XR(3),Xn3l,VR(3),VII3l, IROKI3,2> 


744 


DOUBLE PRECISION fiCOTR , RCCT I , RCCTR ; , RGCT I E , TEMP , TEMP2 ,.\MA x ,C 1 ,C2 

r 


1 

Shl.M.XR.XI.VR.V I.R, ZEPC .CCERR, A 


745 


DOUBLE PRECISION DABS , CS IGN , CSCrT , UMAX l 


745 


INTEGER COLNI, COUNTE, T2 

r ^ 

747 


101 = 1 


748 


103 = 3 


749 


RCOTR = RCCTRt 

sS* 

750 


RCOTI = RCCTIE 

z s 

751 


N = NE 

■ 

752 


MM = M«M - 1 

s S‘ 

753 


M = N - 1 

*' ^ * 

754 


NPl = N ♦ 1 

j: S 

755 


IVCl = IVC - 1 

b. <. • 

7'56 


IVC2 = IVCl - 1 

Z S * 

757 


COUNT = 1 

r.S 

758 


DO 4CC 1=1, N 


'759 


h(I,n=C.CDC 


760 


xRm=o.cDc 


761 

400 

CONTINUE 


762 


CLIH = 1.C6-4 

E S’ 

763 


IF(RQOTI) 1, £3, 1 

5:5’ 


c 


ES 


c 

COMPLEX EIGENVALUE. 

rS' 
' < 

764 

1 

TEMP = - RCCTR - RCCTR 

•7 S 

765 


ISV = 2 

cS’ 

766 


TEHP2»RC!0TR*RCCTR + RCCT 1 *RCCT I 


767 


JJ = 3CG 

= s 

768 


C0~ 6G6 I = i, N 

r 5* 

769 


IFIT2) 600, 6C3, 6C0 

-S* 

770 



CO 602 J = It N 

■ <■. 

771 


JJ « "jJ ♦ 1 

35* 

772 


IFUJ - 251) £02, £01, 6G1 

c • 

773 

60J, 

JJ = 1 

5' 

774 


READ (t'2) ihlLL.i), LL =; l,25'1) 

” s ' 

775 

602 

BII.J) = A(I,J)»TEMP + HIJJ,1) 

rS 

776 


GO TC 665 

- 

777 

603 

B0~'604 j"= 1, N 

C 5 

778 

604 

BtI,J) = AII,J)»TEMP + BII.J) 

c 3 

779 

60S 

atl,l) = BIl.I) + TEKP2 

T ^ 

780 

“ 606 

At 1,1) = Ali.i) - RCCTR 

:5 

781 


IF<T2 .NE. 0) REhINO T2 


782 


GO TC 7C0 


783 

“■ "607 

IFtlCC) 622. 608. 622 



C 


ES 



C 

MATRIX SINGULAR, 

cS 


iri 



139 


c cs 


784 

622 

1FIIVC2) 623, 625, 623 

ES 

785 

623 

DO 624 LL = 1, N 

SS- 

786 


V.(LL,2J=C,CCC 


767 

624 

XItLL)=C.OCC 


788 


IFtIVCl) 625, 514, 625 

ES’ 

789 

625 

DO 626 LL = 1, ,N 


790 


h(LL,4J=C.CCC ‘ ‘ 


791 

626 

VI(LL)=O.CC" 


792 


GO TC 511 

E3‘ 


C 


ES’ 


C 

MATRIX NOT singular. 

es’ 


c 


cS' 

793 

608 

CO 609 LL = 1, N 

es‘ 

794 


ta(LL,l)=l.CCO 


795 


V.(LL,2)=1.:CC 


796 


«(LL,3)=1.C0C 


797 

609 

h(LL,4)=l,000 


798 

699 

IF(I\)C2) 61C, 612, 610 

ES’ 

799 

61C 

CO 611 I = 1, N 

cC' 

800 


12 = IRCW( 1,21 

es’ 

8C1 


XKI2) = ViH.il*RCCTl 

cS’ 

802 


DC 611 J = 1, N 

ES' 

803 

611 

XKI2) = XK12) ♦ Atl,J}«N(J,2) 

ES’ 

804 


IF(IVCl) 612, 500, .612 

6S' 

805 

612 

DC 613 I = 1, N 

ES' 

806 


vnn = nti,3)»RCCTi 

ES' 

807 


DO 613 J =; 1,N 

ES' 

808 

613 

VUI)'^ Vim * A(J,I|«U(J,4) 

ES’ 

809 


GO TO 499 

CS' 

810 

615 

CERR = O.C 

ES' 

all 


CCERR=-J.GDO 


812 


IF(IVCZ) 616, 619, 616 

ES' 

813 

616 

CG 618 I = 1. N 

ES' 

814 


XRU) * -V<<I,2) 

ES' 

815 


CO 617 J =5 1, N 

tS' 

816 

617 

XR(I)-= XR(I) * A(I,JI*XI(J) 

ES’ 

817 

618 

'xR(i) = XRID/RGCTI 

cS' 

818 


IF(IVCl) 61S, 633, 619 

ES' 

819 

619 

CD 621 I = 1, N 

£S' 

820 


VRd) = -H(I,4) 

ES' 

821 


OC 620 J 1, N 

ES' 

822 

620 

VRd) = VRd) ♦ AlJ^r)«VIU) 

PS' 

823 

621 

VR( I ) = VRd l/RGCTI 

ES' 


C 


£ S' 


C 

SEARCH VECrCRS FOR LARGEST ELEMENT ANC NCRKALIZS. 

3S' 


c 


£ S ' 

824 

627 

AHAX>0.000 


825 


DO 629 L = 1, N 

£-S' 

826 


TEMP = VR{L)**2 + VI(L)*«2 

f S' 

827 


IFITEKP - AMAX) 62S, 629, 628 

is' 

828 

628 

AHAX = TEMP 


829 


12 = L 

£S' 

830 

629 

CONTINUE 

ES' 

831 


Cl = VR(I2)/AKAX 

tS' 



1J|0 


832 

833 

834 

835 

836 

630 

837 

838 

839 

631 

840 

632 

841 

633 

042 

843 

844 

845 

634 

846 

847 

635 

848 

849 

850 

851 

852 

853 

636 

854 

855 

856 

637 

857 

C 

C 

C 

638 

85B 

859 

860 
861 
862 
863 

639 

864 

647 

865 

866 

64 2 

867 

640 

868 

869 

641 

870 

871 

644 

872 

373 

64 5 

874 

875 

646 

87o 

877 

873 

64 8 

879 
8 89 
881 

667 

882 

649 

883 



C2 = -VI( I2)/AI^AX 
CG 630 L = 1, .N 
temp - VI (L) 

VIIL) = VP.(L)*C2 + TEPP^Cl 
VRIL) = VR(LI«C1 - TEHP»C2 - 
IFICCUNT .EC. U GC TC 632 
CC 631 LL = I, N 

CCERR=DPAX1 (CCERBfDAPS IVR (LLI-h (LL,3) I , CA P$ t V It LL t-« ( LL , A ) ) ) 

IF{ IVC2I 6 3.3, 636 , 633 

AMAX='J.00C 

CO 635 L = 1, N 

TEKP = XR(LI««2 + XIlLt««2 

IFITEMP - APAXI 635, 635, 63A 

Amax = TEMP 

12 = L 

CONTINUE 

Cl = XR( 121/AMAX ' 

C2 = -XK I2I/AHAX 
00 636 L = I, N 
TEMP = XI (LI 

XKLI = XR(L)»C2 ♦ TEMP*C1 
XR(L) = XR(U»C1 - TEMP*C2 
IF(CCUNT .EQ. 11 GC TC 6A6 
CO 637 LL = I, N 

DCERR=DMAX1(DCERR,DA0S tXR 1LU-HILL,11 1 , C A PS t X I ( LL )-w ( LL , 2 )) 1 

TEST FOR CCNVERGENCE. 

TFICCU.iT ,=C. 1) GC TC 6A6 
CERR=CCERR 

IFICERR .GE. 1.3E-A1 GC TC 639 
IFICERR .GE. CLIP) GC TC 648 
CLIH = CFRR 

IF(CLIR .LE. l.OE-ei GC TC 640 
TFICCLNT .GE. IS) GC TC 68 
COUNT = COLNT 4 1 
IF(RCOTI) 642, 673, 642 
IF(1VC2) 64-:, 644, 64.' 

DO 641 LL = 1, N 
V»tLL,l) = XR(lL) 

WILL, 2) = XT(LL) 

IF( IVCn 644, 610, 644 
CC 645 LL * 1, N 
h(LL,3) = VR(LL) 
h(LL,4) = VI(LL) 

Gr TT 699 
CERR = 0.3 
DC6RR = J.CC.; 

IF(ICC) 643, 647, 648 
ERR = CcRR 
CntNTE = CCLM 
IF{RCGTI) 667, 668, 667 
or 649 I = 1, N 
All, I) = All, I) 4 RCCTR 
RETURN 



i S ' 
OS' 
JS' 
c > 



- c ■ 


ES' 



2 S ' 
OS' 
; c ■ 


£ S ' 



c S ' 

r ' ' 



fcS 


■- ^ 
• t: 


•: c. 




- J 
c 


rS 
- 5 
1 1 


M iA tr tA lA 



1^11 


88A 68 PRINT ICIt RCCTRi RCCTI, CERR tS' 

885 GO TC 6A8 =S' 

C tS- 

C REAL EIGENVECTORS. £.5' 

■ c ■ ■ - - gj. 

886 60 rsh * 1 es’ 

887 _ DC 651 I =>.lj„N . cS' 

888 CO 65C j =5 i, N ES' 

889 650 BII,J) = AU.J) tS' 

8 90_ 6 51 BIJjlIJ = 8(1,1) - RCCTR . . . tS' 

e'91 GO TO 700 ES' 

892 652 IF(ICC) 680, 685, 680 ES' 

C . ES> 

C SINGULAR HATRIX. ES' 

C ES' 

893 680 IFnVC21 681. 683. 68.1. _ ES' 

89A 681 do 682 L =• 1, N ES' 

895 682 XKD^O.OOO 

896 IF(IVCl) 6 83, 51 A. 683 E5' 

897 683 CO 6EA L - 1, N ES' 

898 6EA -71(11=0.000 

899 GO TC 511 . . _ ... E.S' 

C . . ES' . 

C MATRIX NOT SINGULAR. ES' 

C _ _ _ . ES' 

goo' 685 i’F(!VC2) 653, 656, 653 ES' 

901 653 DO 6SA L - 1, N ES' 

902 65A XI(L1 = 1.£D.Q 

903 iF(IVCl) 656, 500, 656 ES' 

9CA 656 DO 657 L = I, N ES' 

905 657 V I (L 1=1,000 „ .. 

906 GO TO AS9 EE' 

C EE' 

C NCRHALIZE REAL VECTORS. _ .ES.' 

■ " "C ' “ ■ ‘ “ ES' 

907 655 CERR = C.C ES' 

908 DCERR=O.Ofi!) _ ... . 

g‘09 IFMVC2) 658, 662, 558 ES' 

910 658 C1=C.CDC 

911 C2=C.gp_0 _ 

912 ' CO 660 L 9 1, N ES' 

913 TeMP=DABStXHL» ) 

91A IF (TEMP - Cl) 660j^660j, .659. . .. . ...EJ' 

915 659 Cl = TEMP ES' 

916 C2 = XKL) ES' 

917 660 CCNTJNUE ES' 

918 *' C'r'661 L ^ 1, N ES' 

919 XKU = Xia)/C2 ES' 

920 OCERR=DMAXl(DCERR,DABS(Xl(U-XR|Lni__ .... 

921 " 661 XR(L) = XI(L) ES' 

922 IF(lVCl) 662, 63£, 662 ES' 

923 662 C2 = O.ODD _ . .. . 


925 CO 66A L = I, N ES' 

926 . TEHP = DA8S(.VI (L)J__ 



. 1.42 


927 


IF(TEMP - Cl) 664, 664, 663 

tS> 

928 

663 

Cl = TEMP 


929 


C2 = VI {L) 

£S' 

930 

664 

CONTINUE . . . 

es' 

931 


cn 665 LL = 1, N 

es> 

932 


VKLU = VI(tL)/C2 

ES> 

933 


DCERR*CMAXl(nCEllR,CAeStVI ILL)-H(LL,U U 


934 


h (i.L,l)=VHLL ) 


935 

665 

VR(LL)=MLL,1) 


936 


GO TC 638 ' . ■ _ 


937 

668 

iF(IVC2) 669, 671, 669 

ES^ 

938 

669 

00 670 L » 1, N 

ES' 

939 

670 

XI(U«O.OOC 


940 


IF(IVCl) 671, 70, 671 

c c> 

941 

671 

CO 672 L = 1, N 

is' 

942 

672 

VI(L)=O.CDC 


943 

70 

RETURN 

ES' 

944 

673 

IF(IVC2) 674, 502, 674 

ES' 

945 

674 

CO 675 I =« 1, N 

■ es.' 

946 


12 = IRCV,( 1,2) 

ES' 

947 

675 

XKI2) = XR(l) 

uS' 

948 


GO TC 5CC 

cS ' 


C 


ES' 


C 

BACK SURSTITUTTCN SECTION, 

ES' 


C 


... . ES' 

949 

499 

IF(IVC2) 5CC, 502, 5CC 

ES' 

950 

500 

DO 501 1 = 2, N 

cS' 

95X 


11 = I - 1 

is' 

952 


DO 5C1 J = 1, 11 

ES' 

953 

501 

xim = xim - E<i,j)»xiu) 

ES' 

954 

511 

IFUVCl) 502, 514, 502 

tS' 

955 

502 

00 510 I = 1, N 

cS- 

956 


11 = I - 1 

ES' 

957 


IFlIl) SC3, 5v5, 5C3 

cS' 

958 

503 

CO 504 J = 1, 11 

ES’ 

959 

50.4 

vim = Vim - eu,i)*vi(j) 

■ES' 

960 


IFUCC) 505, 506, 505 

ES' 

961 

5C5 

IHB(I,in SC6, 5C7, 5C6 

ES’ 

962 

5C6 

Vim = vim/Bn,i) 

ES' 

963 


GC TC 51G 

c-S' 

964 

507 

IF(Vim) 508, 5C9, 5C8 

ES' 

965 

508 

Vim = vim«i.0E+i5 

ES' 

966 


GO TO 510 

ES' 

967 

509 

viii) = i.c 

£S' 

968 

51C 

CONTINUE 

lS' 

969 


IF(IVC2) 514, 525, 514 

cS’ 

970 

514 

DC 522 I = 1, N 

ES' 

971 


IR = NPl - I 

ES' 

972 


IFII - 1) 515, 517, 515 

ES’ 

973 

515 

12 = IR 4 1 

ES' 

974 


cn 516 J = 12, N 

ES' 

975 

516 

XI(IR) = Xltlft) - BIIR,J)4XI(J) 

ES' 

976 


IF(ICC) 517, 518, 517 

ES' 

977 

517 

IFIBdR.IR) ) 518, 519, 518 

cS' 

978 

518 

XIIIR) = XKIRI/BIIR.IR) 

ES' 



1^3 


979 


60 TC 522 

es' 

980 

519 

tFCXIlIRl ) S2C, 521, 520 

E'^ 

981 

52C 

XntR) = XmR)*l,CE+15 

ES' 

902 


.SC TC.522_ 



983 

S21 

XKIRl^l.ODQ 


98<i 

522 

CONTINUE 

ES' 

985 


IFdVCn 525, S29, 525 

£S' 

986 

525 

DO 526 I » 2, N 

£S' 

987 


IR = NPl - I 

tS^ 

988 


12 = IR + 1 

JE.S' 

989 


DO 526 J 9 12, N 

ES' 

990 

526 

VKIRl > VMIR) - BU,IR)«VI(JI 

ES' 

991 


C0..52V.U - . 

.... . .... .ES' 

992 


12 = IRGHIL,!) 

E5‘ 

993 

527 

VR(I2) = VI (L) 

ES' 

994 


CO 528 L * 1. N 

ES' 

995 

528 

VKL) - VftlL) 

ES' 

996 

529 

IF(RCOTI) 615, 656, 615 

ES' 


C 


ES' 


C 

FACTOR FATRIX. 

ES' 


C 


ES' 

997 

700 

ICC = 0 

ES' 

998 


Shl»l,0D72 ^ 


999 


DO 701 LL « 1, M 

ES' 

1000 

701 

IROWILL.II ■ LL 

ES' 

1001 


DO 70S K ^ 1, N1 

ES' 

1002 


AHAX*0ABS{8IK,K> ) 


1003 


IKAX * K 

ES' 

1004 


K1 » K ♦ 1 

G5' 

1005 


DO 702 I 5 XI, N 

ES' 

1CC6„. 


IFlAHAX.GT.CABSIBfl.Xn > GO TG 702 


1007 


AHAX»DA8S(a(t,KI > 


1008 


IMAX = I 

ES' 

1009 

702 

CONTINUE 

ES' 

1010 


IFtAXAX .LT. SUL) SMI » AFAX 

ES' 

1011 


IF(AFAX.GE*1.0D-25I GO TC 723 


.1012 


BIK.KMO.OOO 

• 

1013 


ICC = ICC ♦ 1 

ES' 

1014 


GO TO 7C8 

ES' 

1015 

723 

IFI FFAX .EC. Kl GO TC 704 

ES' 

1016 


DO 703 J ^ 1, N 

ES' 

1017 


ANAX • 8IR*J) 

ES' 

1018 


BlX.J) - 

ES' 

1019 

703 

B(IMAX,J) > AFAX 

ES' 

1020 


12 > IRQMIK.I) 

ES' 

1021 


IROWCK,!) > IROWIIFAX.l) 

ES' 

1022 


IRQH(IM*X,U > 12 

ES' 

1023 

704 

DO 707 t 4 XI, N' 

ES' 

1024 


iFfBIUKM 705 . 707 . 705 

ES' 

1025 

705 

B(I,K) < B(I,K)/B(K,X) 

ES' 

1026 


DO 706 J - XI, N 

ES' 

1027 

706 

B(1,J) > B(I.U) > BIX.JIfBII.KI 

ES' 

1028 

707 

CONTI NUi 


1029 

708 

CDMTINU2 ' 


1030 


ANAX-0A«2»lH»»im ^ ^ 




i44 


1C31 IF(AKAX-1.0D-25) 712,712,713 

1032 712 B(N,N)=0.000 

1033 SHl=C.ODO 

1034 ICC » ICC ♦ 1 E 

1035 GO fC tog 6 

1036 713 IF(APAX .LT. SHI) SUl > AMX E 

1037 709_iJF < liX . LEj,. ISH) GC TC UQ £ 

1038 IFlMr) 1050,1050,1051 

1039 1051 HRITEI 103,102) ICC 

1040 OUN 1E_5_C ' • c 

1041 RETORN E 

1042. 1050 HR1TE(I03,1QS2) ICC 

1043 710 00 tJXLL « 1, N _ £ 

1044 12 X IRGH TlL,1) c 

1045 711 IRCH(I2,2) * LL 6 

1046 I FJ RCDTI) 6C7, 652, 607 ^ „ - 


runmmiif/r/ian hakriki* ,• :>udkuui inc tibVbi. t 

IFOUNC AN EIGENVALUE CE APPARENT PULTIPL ICITY* , E 


, CCHPUTATION OF EIE 


1048 

■iyAg" 

1050 


2GENVECTQR(S ) CONTINUES AT USER S OPTION*//) 

101 FORHATt 38H0P0RE THAN 15 LCCPS FOR EIGENVECTOR OF,2E12-4, 
.2 14H difference CF.E12. 4) . . . 


102 FORMATa6HO«**«MAt»AIHC»**»' • 14. 71H ZERCS ON DIAGONAL Of FACTOREDc 
1 HATRIX* CHECK FOR PULTIPLE EIGENVALUES. /20X, E 

2» SUBROUTINE EIGVtC HILL NOT P ERFORM CONPUTATION FCR THIS EIGENVECE 
3T0R *//> E 

END 


/DATA 


</' M f/\ »/4 tfi [/» 1/1 I/V tA (n I/I 1/1 t/l 4/1 f/l </) tA tA 



